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ABSTRACT 


During  the  past  two  years  (12/00-3/03),  we  have  developed  a  very  efficient  and  reliable 
direct  numerical  simulation  (DNS)  and  large  eddy  simulation  (LES)  software  (DNSUTA) 
for  compressible  flow  in  a  curvilinear  coordinates  which  is  a  continuation  of  our  previous 
effort  supported  by  AFOSR.  The  software  includes  a  high-quality  grid  generation  code 
which  can  generate  smooth  and  near-orthogonal,  body-fitted  grids  and  a  parallel  efficient 
flow  solver  which  applies  sixth-order  compact  scheme,  eighth-order  filter,  LU-SGS  flow 
solver,  non-reflecting  boundary  condition,  high-order  weighted  compact  scheme  for 
shock  capturing,  structure  function  sub-grid  model,  and  MPI  parallel  computation.  The 
code  has  been  successfully  used  for  prediction  of  flow  transition  around  flat  plate  and 
airfoils  with  various  attack  angles.  The  code  has  also  been  used  for  flow  separation 
control  for  low  Reynolds  airfoils  with  steady  and  pulsed  blowing  jets.  All  of  these 
simulations  are  real  time  dependent,  3-D,  with  general  geometry  and  have  demonstrated 
qualitative  or  quantitative  agreement  with  experiment  or  theoretical  work.  We  believe  the 
topics  described  above  are  very  challenging  and  our  accomplishment  is  very  significant 
for  both  basic  science  and  Air  Force  mission.  The  high-order  weighted  compact  scheme 
and  non-reflecting  boundary  conditions  in  a  curvilinear  coordinate  are  the  original  work 
conducted  by  us  which  is  a  significant  contribution  to  the  CFD  community. 

The  code  is  currently  used  by  NASA  Langley  for  prediction  of  flow  transition  and 
turbulence  through  a  technology  transfer  contract.  The  code  has  passed  all  NASA  test 
cases  and  demonstrated  its  high  accuracy  and  high  reliability.  The  code  and  related  work 
have  been  recognized  by  NASA  Langley  Researchers  as  unique  in  the  United  States.  The 
code  also  attracts  significant  interests  from  US  Navy  for  study  of  wakes  and  propeller 
flows. 

Because  of  the  unexpected  notice  issued  by  AFOSR  that  the  Optional  2  of  our  grant  will 
not  be  exercised  due  to  funding  restrictions,  we  lost  the  optional  fund  for  the  third  year 
and,  therefore,  we  did  not  have  enough  time  to  complete  all  of  work  we  proposed. 


3 


Chapter  1 
Introduction 

Turbulent  flow  prediction  and  analysis  are  very  important  to  Air  Force  missions.  There 
are,  in  general,  three  ways  to  conduct  turbulent  flow  prediction  -  direct  numerical 
simulation  (DNS),  large  eddy  simulation  (LES),  and  Reynolds  averaged  Navier-Stokes 
solver  (RANS),  while  the  time-dependent  Navier-Stokes  equations  are  still  widely 
accepted  as  a  governing  system  for  both  laminar  and  turbulent  flows. 

The  classical  averaging  method  is  so-called  Reynolds  averaging  which  appeared  around 
one  hundred  years  ago  and  is  still  widely  used  by  industry.  RANS  is  an  economic  way  to 
simulate  turbulent  flow  assuming  a  mean  flow  followed  by  a  perturbation.  RANS  only 
studies  the  flow  mean  value.  However,  in  many  important  flow-related  topics  including 
boundary  layer  receptivity,  flow  transition,  flow  separation,  wakes,  interaction  of  flow 
separation  and  wakes,  noise  generation,  vortex  dynamics,  flow-structure  interaction, 
unsteady  turbulence,  turbulent  heat  transfer,  hypersonic  chemically  reacting  flow,  shock- 
turbulent  boundary  layer  interaction,  suppression  of  turbulence,  etc.,  which  Air  Force  is 
highly  interested  in,  only  the  instantaneous  quantities  are  critical.  For  these  Air  Force 
interests,  DNS  or,  at  least,  LES  must  be  used.  Although  DNS  is  very  expensive,  there  is 
still  a  large  demand  for  developing  efficient  DNS  technology. 

The  RANS  system  is  not  self-closed  and  many  correlation  terms  need  to  be  modeled. 
Most  popular  eddy  viscosity  models  assume  Reynolds  stress  is  proportional  to  strain  of 
the  mean  flow.  Unfortunately,  modem  fluid  mechanics  study  has  found  there  is  no  such  a 
correlation.  The  interaction  between  large  length  scales  and  small  length  scales  is  very 
complicated.  Most  early  turbulence  models  assume  the  turbulence  is  isotropic,  an 
assumption  that  is  now  understood  to  be  unrealistic  especially  for  wall-bounded  flow. 
Some  new  turbulence  models,  like  Reynolds  stress  models,  have  been  offered  in  attempts 
to  remove  the  assumption,  but  these  models  introduce  higher-order  closure  problems. 
RANS  is  currently  used  widely  by  industry  and  is  expected  to  be  relied  upon  until  we 
have  new  reliable  and  feasible  approaches.  However,  RANS  has  fundamental  limitations 
and  deficiencies.  Turbulence  models  are  usually  empirical  and  case-related.  In  general, 
they  are  good  for  some  cases  largely  because  of  adjustment  of  coefficients,  but  not 
accurate  for  other  cases.  It  is  really  hard  or  almost  impossible  to  find  a  universal  model 
which  can  predict  well  for  all  flow  cases.  As  pointed  out  by  AGARD-CP-551  at  its 
Foreword,  “one  of  the  main  conclusions  was  that  no  one  model,  or  type  of  model,  could 
predict  all  the  test  cases  to  good  engineering  accuracy.  The  implication  of  this  result  is 
that  for  certain  cases  Large  Eddy  Simulation  or  Direct  Numerical  Simulation  may  be  the 
only  recourse  for  obtaining  an  accurate  prediction”  (AGARD,  1994). 

Then,  we  have  to  look  back  to  the  original  time-dependent  Navier-Stokes  equations  in 
our  search  for  solutions.  This  leads  us  to  direct  numerical  simulation  (DNS).  DNS  is  the 
only  numerical  approach  dealing  with  the  exact  solution  of  3-D  time-dependent  N-S 
equations  without  any  ad  hoc  assumption  and  in  which  leads  all  cases  to  good 
engineering  accuracy.  However,  it  is  still  too  expensive  if  we  want  to  resolve  all  length 
scales  including  the  Kolmogorov  scale  especially  for  high  Reynolds  number  and  general 
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geometries.  We  need  to  refine  the  grids  in  certain  areas  where  the  small  length  scale  is 
important.  More  realistic  at  present  time  is  to  look  for  some  compromise  which  resolves 
most  significant  large  vortex  and  leave  small  length  scales  (smaller  than  the  grid  size)  for 
modeling.  This  will  lead  to  large  eddy  simulation  (LES).  In  LES,  the  time-  and  space- 
dependent  large  scales  will  be  simulated  and  models  are  only  needed  for  scales  smaller 
than  the  grid  size.  Although  the  small  length  scales  are  more  dissipative  and  easier  to 
model,  the  grid  size  must  be  relatively  small  until  the  subgird  scale  is  less  important  in 
LES  in  order  to  achieve  accurate  results. 

DNS  for  flows  with  low  Reynolds  number  and  simple  geometry  has  been  successfully 
performed  and  it  is  widely  recognized  that  DNS  can  help  understand  flow  physics  and 
can  check  and  improve  turbulence  models.  The  question  is  whether  or  not,  in  any  sense, 
DNS  can  work  for  real  world.  The  fundamental  problem  with  DNS/LES  seems  to  be 
solved  and  our  experience  shows  DNS  can  be  used  for  flow  transition  and  LES  can  be 
used  for  fully  developed  turbulent  flow  for  a  spectrum  of  engineering  applications. 

The  recent  developments  in  computer  performance  and  numerical  algorithms  seem  to 
bring  hope  to  the  DNS/LES  community  that  DNS/LES  can  in  the  near  future  be  modified 
for  high  Reynolds  numbers  and  complex  geometries.  There  are  several  hopeful  signs: 
First,  computer  capability  has  increased  dramatically  every  year.  The  new  generation  of 
parallel  machines  with  1000  processors  and  distributed  memory  can  have  Terra-bytes  of 
memory.  The  price  of  memory  has  dropped  sharply.  The  problem  with  computer  memory 
requirements  seems  to  be  gradually  disappearing.  The  CPU  capacity  remains  a  challenge. 
The  fastest  parallel  machine  with  1000  processors  nowadays  can  perform  over  two  Terra- 
flop,  raising  by  one  order  the  Reynolds  numbers  that  can  be  covered.  Nowadays,  at  least 
200  millions  of  grids  can  be  used  for  flow  simulation  on  currently  available 
supercomputers.  More  recently,  Linux  PC-Cluster  has  been  well  developed  and  widely 
used  which  reduces  the  computer  cost  by  one  order  to  the  range  of  $20-30k  for  16  CPUs, 
which  is  affordable  to  most  aerospace  companies  and  encourages  more  scientists  and 
engineers  to  use  DNS/LES  for  flow  simulations.  Second,  DNS  for  complex  geometry  is 
feasible.  It  requires  a  high-order  grid  generation,  high  grid  quality  and,  moreover,  high 
order  discretization  and  high  order  filter  for  general  curvilinear  coordinates.  A  typical 
turbulent  flow  has  a  transition  process  which  includes  receptivity  of  environmental 
disturbances,  linear  growth,  nonlinear  instability,  breakdowns,  and  transition,  eventually, 
to  a  fully  developed  turbulent  flow.  The  length  scale  changes  at  every  stage.  In  the  main 
flow,  the  length  scale  can  be  thought  of  as  0(1) .  In  general,  the  receptivity  and  early 

transition  processes  are  dominated  by  large  scales  of  the  order  of  0(Re)’2 .  At  these 
stages  the  grid  size  need  not  to  be  very  small.  Probably  it  is  small  enough  if  it  can  resolve 
the  Tollmien-Schlichting  waves.  Therefore,  DNS  can  be  used  for  main  flow  and  early 
transition  stages  for  some  realistic  Reynolds  numbers  and  realistic  geometries.  When  the 
flow  enters  the  non-linear  transition  zone  and  finally  becomes  fully  developed  turbulent 
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flow,  the  smallest  length  scale  may  reach  0(Re)“4 .  When  the  Reynolds  number  is  very 
high,  we  still  do  not  have  enough  grids  to  resolve  them.  We  can  only  resolve  those  scales 
greater  than  or  equal  to  the  grid  size  and  have  to  model  other  flows  with  high  Reynolds 
number  and  general  geometries.  DNS/LES  for  complex  geometry  was  criticized,  but  now 
an  increasing  number  of  researchers  (e.g.  Liu,  et  al,  1994,  1998a,  1998b;  Jiang,  1999a, 
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2001b;  Shan  2000a,  2000b;  Visbal  et  al.,  1998;  Adams,  2000;  Rizzettz  et  al,  2001;  etc.) 
are  working  on  this  direction.  Nowadays,  DNS/LES  can  capture  the  major  feature  of 
transitional  and  turbulent  flow  around  a  3-D  airfoil  based  on  the  currently  available 
parallel  computers. 

Looking  around  for  the  trend  of  CFD,  we  can  find  more  and  more  scientists  have  realized 
the  problem  with  the  Reynolds  averaging  and  focused  their  research  on  DNS  or  LES,  The 
time  has  come  to  orient  numerical  simulation  to  the  real  world.  The  “real  world”  for 
flows  around  flight  vehicles  has  the  following  features: 

(1)  The  flow  is  real-time  dependent, 

(2)  The  vehicles  are  geometrically  very  complex, 

(3)  The  vehicles  operate  at  Mach  numbers  which  range  from  high  subsonic  through  high 
supersonic,  or  even  hypersonic  in  the  near  future,  where  finite-rate  chemical  reactions 
must  be  considered, 

(4)  The  flow  field  is  a  combination  of  laminar  flow,  transitional  flow,  and  fully  developed 
turbulent  flow. 

It  is  understood  that  the  following  challenges  for  DNS/LES  of  turbulent  flows  should  be 
met: 

(1)  High-order  accuracy  discretization  in  both  time  and  space  in  a  curvilinear  coordinate, 

(2)  High-order  shock-capturing  for  transonic  and  supersonic  flows. 

(3)  High-order  grid  generation  and  grid  mapping  for  general  configurations. 

(4)  Very  fine  resolution  for  boundary  layers,  especially  the  near-wall  region,  to  resolve 
small-scale  eddies  involved  in  the  transition  process  and  turbulence. 

(5)  High-order  filters  and  reliable  subgrid  scale  models. 

(6)  High  code  efficiency  to  reduce  CPU  cost  to  a  reasonable  level. 

These  are  very  challenging  requirements.  However,  our  approach  has,  in  principle, 
provided  a  solid  base  to  meet  the  above  requirements  and  is  therefore  a  feasible  approach 
to  DNS/LES  for  engineering  applications. 

The  serious  problem  is  that  the  small  length  scales  will  become  important  when  the  flow 
enters  the  late  transition  stage  (non-linear  and  breakdown  stages).  In  the  case  that  the 
Reynolds  number  is  high  and  geometry  is  complicated  where  our  grid  resource  still 
cannot  provide  enough  resolution,  we  can  then  use  LES.  Therefore,  we  can  say  that 
based  on  the  capability  of  current  available  supercomputers,  we  can  use  DNS/LES  to 

solve  turbulent  flow  in  turbomachinery  which  has  a  Reynolds  number  of  around  105  . 
The  DNS  can  be  also  used  for  3-D  airfoil  flow  transition  for  a  Reynolds  number  of 

around  105  -106 .  The  Reynolds  number  can  be  near  one  order  higher  if  we  use  subgrid 
models  (LES). 

Our  work  has  considered  both  the  basic  science  advances  and  the  Air  Force  needs.  It  is 
therefore  anticipated  that  the  success  of  our  NDS/LES  code  will  provide  the  following 
potential  to  AFOSR  and  other  federal  agencies  (NASA,  ONR,  etc): 
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(1)  Provide  a  detailed  description  of  the  whole  process  of  boundary  layer  receptivity, 
flow  separation,  wakes,  interaction  of  flow  separation  and  wakes,  noise  generation, 
vortex  dynamics,  shock-turbulent  boundary  interaction,  suppression  of  turbulence,  flow 
transition  and  turbulence  structure.  This  will  have  a  notable  impact  on  fundamental  fluid 
mechanics. 

(2)  Help  understand  the  origin  of  turbulence  such  as  how  the  environmental  disturbance 
in  free  stream  induces  unstable  wave  modes  in  boundary  layers  and  what  causes  the  non¬ 
linear  instability  and  flow  breakdown  and  how  the  unstable  wakes  induce  the  original 
inviscid  shear  layer  instability  for  airfoils  with  attack  angle  (see  Jiang  et  al,  2003). 

(3)  Provide  a  “Computational  Wind  Tunnel”  which  can  give  a  full  simulation  of  turbulent 
flow  around  airfoils  for  flow  transition,  separation,  and  flow  control. 

(4)  Provide  assistance  to  transition  control,  drag  reduction,  noise  reduction,  and 
improvement  of  aircraft  design.  As  advances  are  made  in  supercomputer  performance 
DNS/LES  for  engineering  design  will  become  feasible  with  more  complex  geometry  and 
high-Reynolds  number.  It  will  provide  a  right  answer  to  flow  prediction  on  drag  and  heat 
transfer. 

During  the  past  two  years,  12/00-3/03,  we  have  successfully  developed  a  very  accurate 
and  reliable  DNS/LES  code  (DBNSUTA).  The  code  has  been  successfully  used  to 
simulate  the  whole  process  of  flow  transition  for  flat  plate  and  airfoils  at  both  low  speed 
(Af„  =0.2  for  flat  plate  and  airfoils)  and  high  speed  (flat  plate  at  =4.5).  The  DNS 
and  LES  results  have  been  well  documented  and  validated  (see  Liu  et  al,  1991a,  1991b, 
1993,  1994a,  1994b,  1995,  1998a,  1998b;  Shan  et  al,  1999a,  1999b,  1999c,  2000a,  200b; 
Jiang  et  al,  199a,  199b,  2001a,  2001b,  2003;). 

The  code  has  been  validated  by  the  authors  in  a  number  of  publications  and  by  NASA 
users.  NASA  Langley  personnel  (Choudahari  et  al,  2003a  and  2003b)  used  DNSUTA  for 
a  number  of  cases  of  flow  transition  including  flat  plate,  swept  wing  and  straight  cone 
and  found  the  validation  is  very  satisfactory.  The  code  also  attracts  significant  interests  to 
US  Navy  for  study  of  wakes  and  propeller  flows.  As  a  Navy  David  Taylor  researcher 
points  out  (personal  communication),  we  have  no  sense  to  wait  for  5  years  for  developing 
LES  code  for  general  geometry,  but  not  to  use  ‘DNSUTA’  to  simulate  wakes  now. 

Recently,  Dr.  Liu  and  his  colleague  have  use  DNSUTA  to  successfully  simulate  the  flow 
separation  and  transition  around  airfoil  NACA0012  with  four  degree  of  attack  angle  and 
Reynolds  number  of  100,000.  They  also  use  DNSUTA  for  flow  control  with  pulsed  and 
screwed  blowing  jets.  These  achievements  clearly  show  the  purpose  of  this  project  has 
been  achieved  successfully.  Unfortunately,  we  cannot  complete  all  tasks  we  proposed  due 
to  the  sudden  stop  of  Option  2  fund  noticed  by  AFOSR  contracting  office.  We  lost 
Optional  2  fund  for  the  third  year  and  then  do  not  have  enough  time  to  complete  all  work 
we  proposed. 
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Chapter  2 


Our  Technical  Approaches 

Our  DNS/LES  code  adopts  state  of  art  approaches  including  orthogonal  grid  generation, 
high-order  Jacobien  in  curvilinear  coordinates,  high-order  compact  scheme,  high-order 
filter  (filtering  the  LES  solution  at  each  time  step),  high-order  weighted  compact  scheme, 
non-reflecting  boundary  condition,  efficient  Navier-Stokes  solver,  effective  dynamic  sub¬ 
grid  model,  and  parallel  computation.  The  high-order  weighted  compact  scheme  is 
particularly  critical,  which  can  capture  the  shock  very  accurately  and  keep  high-order  for 
smooth  areas.  The  effectiveness  of  the  subgrid  model  is  another  critical  factor  that  would 
ensure  the  accuracy  of  LES  with  rather  coarse  grids.  Our  work  on  weighted  compact 
scheme  and  non-reflecting  boundary  condition  in  a  curvilinear  coordinate  are  original  and 
are  a  significant  contribution  to  the  CFD  community. 

2.1  Contravariant  velocity-based  governing  equations: 

The  three-dimensional  time-dependent  Navier-Stokes  equations  in  generalized 
curvilinear  coordinates  (£,  tj,C)  can  be  written  in  conservative  forms: 

1  B<2  3{E-E,)  3{F-Fr)  d{F-F,) 

J  1  drj  d( 

The  vector  of  conserved  quantities  Q ,  inviscid  flux  vector  (E,  F,G),  and  viscous  flux 
vector  (Ev ,  Fv ,  Gv )  are  defined  via 


' P  N 

(pu  1 

(pV 

' pw  N 

pu 

’E  =  l 

pUu  +  p£x 

■f-7 

pVu  +  pTJx 

,  G  =  — 
J 

pWu  +  pCx 

Q  = 

pv 

pw 

pUv  +  p%y 
pUw  +  p£z 

PVV+pTJy 
pVw+  pT)z 

pWV+  pCy 
pWw+  p£z 

Ke  j 

KU(e  +  p )  j 

KV  (e  +  p)  J 

KW(e  +  p )  , 
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'o  N 

„  1 
E^1 

TJx  +Tyxty+rJz 
Txytx  +Tjy  +*Jz 

f’=7 

rxxnx+ryxr]y+rant 

rxyTJx+TyyT}y+TzyTJz 

*Jx  +*Jy  +  TJz 

TxzVx  +*yz‘ny  +?uVz 

Kqx^x+qy4y+Qz^z  J 

^xVx+qyVy+qz1 Iz  , 

Tier Cx  ^  yxC  y  't'zxC z 

GV=J  TxyCx+TyyCy+^Cz 

I'xzCx  "*"  ^  yz  ^  y  ^ zz^>  z 

^IxCx  ^lyCy  GZCZ  j 

where  J  =  js  Jacobian  of  the  coordinate  transformation  between  the  curvilinear 

d{x,y,z) 

(£,  rj,  C)  and  Cartesian  ( x,y,z )  frames,  and  %x,%y,%z,Vx,Vy,Vz>Cx’Cy’Cz  are  coordinate 
transformation  metrics.  The  contravariant  velocity  components  U,V,W  are  defined  as 
U  sk£  +v£  +w£,  V  =  UT]X  +vrjy+wT}z,  W  =uCx  +vCy  +wCz, 

e  =  —  +—p(u 2  +v2  +  w2). 
y- 1  2 

The  components  of  the  viscous  stress  tensor  are  denoted  by  ,  Tyy  ,Tzz,rxy,  Txz ,  ryz .  The 
density,  velocity  components  and  pressure  are  defined  by  p ,  u,  v,  w  and  p. 

2.2  Orthogonal  grid  generation  and  high  order  Jacobian 

In  order  to  achieve  high-order  accuracy,  say  fourth-order  or  higher,  first,  a  high-order 
description  of  Jacobian  metrics,  e.g.,  £x,£y etc.  is  required.  Second,  the  high  quality  of 

grids,  i.e.  continuity  in  curvature  must  be  guaranteed.  Note  that  there  are  several  Jacobian 
grid  mapping  terms,  J,  ...  in  our  governing  equations  which  need  to  be  evaluated 

by  a  high-order  finite  difference  scheme  to  guarantee  the  accuracy  of  our  high-order 
discretization.  As  shown  by  numerical  experiment,  any  discontinuity  or  disturbance  in 
surface  curvature  will  result  in  a  totally  different  output  with  same  inflow  disturbance. 
Accurate  LES  requires  a  very  smooth  gird  generation  where  the  second-order  derivative 
is  continuous  and  the  third-order  derivative  exists  everywhere.  Also  an  orthogonal  grid 
especially  near  the  solid  surface  is  preferred.  A  high-quality  elliptic  grid  generation  has 
been  carried  out  for  2-D  airfoil  and  3-D  Delta  wing  (Shan  et  al,  2000,  Figure  1  and  2).  An 
elliptic  grid  generation  method,  first  proposed  by  Spekreijse  (1995),  is  used  to  construct 
the  three-dimensional  grids.  This  method  is  based  on  a  composite  mapping,  which 
consists  of  a  nonlinear  transfinite  algebraic  transformation  and  an  elliptic  transformation. 
The  grids  are  orthogonal  on  the  delta  wing  surface.  The  three-dimensional  grid  sketch  can 
be  found  in  Figure  2  which  shows  we  have  generated  smooth  and  near  orthogonal  high- 
quality  grids. 
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Figurel.  High-quality  grids  for  2-D  airfoil  Figure  2.  High-quality  grids  for  3-D  Delta 

wing 

2.3  High-order  compact  scheme  and  high-order  filter 

The  governing  equations  are  discretized  in  time  by  implicit  method  based  on  second 
order  Euler  Backward  scheme.  The  sixth-order  centered  compact  difference  scheme  is 
used  for  spatial  derivatives.  High-order  compact  filter  is  employed  to  reduce  numerical 
oscillation. 


In  Eq.  (2.1),  a  second  order  Euler  Backward  scheme  is  used  for  time  derivatives,  and  the 
fully  implicit  form  of  the  discretized  equations  is  given  by 

3Qn+1  -  4  Qn  +  Q"~l  d(En+]  -  £;+1 )  a(F"+'  -  Fvn+1 )  d(cn+l  -  Gf  )_n 

2  JAt  a#  drj  dC  {  } 

Qn+I  is  estimated  iteratively  as: 

Qn+l  =  Qp  +  SQP 
where, 

SQP  =Qp+'-Qp 

At  step  p  =  0 ,  Qp  =  Qn ;  as  5QP  is  driven  to  zero,  Qp  approaches  Qn+l .  The  flux  vectors 
are  linearized  as  follows: 

En+l  ~  Ep  +  APSQP 
Fn+'  ~  Fp  +  BP5QP 
Gn+I  ~GP  +  CP6QP 
So  that  Eq.  (2.2)  can  be  written  as: 

'3 


-1  +  A  tj{D(A  +  D^B  +  D(c) 


R 


(2.3) 


where  R  is  the  residual: 
'3 


R 


1  N 

-Qp -2  Qn  +-Qn~l 
2  2 


■  A tJ [D((E-Ev)+Dn{F-Fv)+D({G- Gv )]p  (2.4) 


Dc,DlVD,  represent  partial  differential  operators,  and  A,  B,  C  are  the  Jacobian  matrices 
of  flux  vectors: 


A-M. 

A“32’ 


B  = 


3F 

dQ’ 


G  =  & 
dQ 


The  right  hand  side  of  Eq.  (2.3)  is  discretized  using  sixth-order  compact  scheme  (Lele, 
1992)  for  spatial  derivatives,  and  the  left  hand  side  of  the  equation  will  be  solved  by  LU- 
SGS  method  (Yoon  et  al,  1992). 


The  general  form  for  compact  finite  difference  schemes  (Lele,  1992)  is  given  by 
fi-f' j-2+a  f' j+a+f' j+i+fi+f' j-2  +a_fj_l+c  fj  +  a+f +  b+fj+2) 


For  sixth-order  compact  scheme: 

P-  =0,  a_  =  pfft  =\*P+  =0’ 


(2.5) 


b  = - ,  a  = — ,  a , 

36  ‘  9 


36 


Because  of  the  use  of  the  central  difference  scheme,  two-point  "saw-tooth"  oscillations 
will  generally  be  generated,  especially  in  the  low  viscous  region.  Those  oscillations  can 
induce  some  spurious  non-physical  waves.  To  avoid  this  phenomenon,  we  need  to  use 
high-order  (8th  or  10th  order)  filter  for  the  solution,  which  is  a  different  concept  from 
filtering  the  equation  in  LES  but  filtering  the  solution,  at  each  time  step.  A  compact  sixth- 
order  filter  can  be  written  as  follows: 


P  fj-2  +  a  fj- 1  +  fj  +  a  fj+l  +  P  fj+2  = 

d 

2  ^  j-j  '  2 


^  +  4  (/,+ 3  +  fj-2  )  +  “(  fj*  2  +  fj-2  )  +  |  (/j+1  +  fj- 1 ) 


(2.6) 


where  /  is  filtered,  f  is  the  original,  and 


a  =  —  fi- 

8, 


3-2 a 
10 


2  +  3  a 


a  =- 


,b  = 


6  +  la 


6+a  ,  2-3 a 

-  d  - - 

20  40 


(2.7) 


2.4  Efficient  Navier-Stokes  Solver 

We  use  LU-SGS  method  (Yoon  et  al,  1992)  to  solve  the  Navier-Stokes  equations. 

The  equation  can  be  written  as 

[^l+&tJ(D;A  +  DvB  +  D(C)]SQp  =R  (2.8) 

The  right-hand-side  of  Eq,(2.8)  is  discretized  using  the  fourth-order  compact  scheme  for 
spatial  derivatives.  The  left-hand-side  of  Eq.(2.8)  is  discretized  using  the  LU-SGS 
method.  The  Jacobian  matrices  of  flux  vectors  given  by  Eq.(2.8)  are  split  as: 


A=A+  +  A~,  B  =  B++B~,  C  =  C++(r 
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where 

= 

l[A±r^7] 

B±  = 

B±tbI] 

= 

\lC±rcI] 

and 

ta  —  K?naa:[|A(j4)|]  +  v 
rB  —  «maz[|A(B)|]  +  v 
rc  =  Kmaz[|  A(C)|]  +  v 


X(A),  X(B),  X(C)  are  eigenvalues  of  A,  B,  C.  The  effects  of  viscous  terms  are  taken  into 
account  by  adding  a  constant  given  by 


r  P  4  t1  i 

v  —  max\- - r—z - ,  —  — 

1  (7  —  l)M?RePr  3  Re  (2  9) 

The  first-order  upwind  finite  difference  scheme  is  used  for  the  split  flux  terms  in  the  left- 
hand-side  of  the  Eq.(2.8) 

[|/  +  A<J(rA  +  rB  +  ro)/|«<3?ji4= 

-aw  [  .4-aor+u,*  - 

+  -  B^SQ^ 

+  CSQl^-CUQ^  ] 


Eq.(2.10)  is  divided  into  the  following  three  steps  and  solved  by  iteration. 
Step  1: 

=  [|/  +  AtJ(rA  +  rB  +  rc)f]_1-R?jfe 

step  2: 

SQij,  k  =  +  [|f  +  AdJ  {rA  +  rB  +  rc)/]-1 

x{AtJ(A+£QUj,k  +  B+6Qlj_hk+C+6Q^k_1)] 

step  3: 
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Vi j,k  =  SQ'.J,l  “  [f 7  +  AtJ (r-t  +  rB  +rc)]  1 

x[AiJ(A-«Q?+liiii  +  BSQl^  +  CSQlM1)} 


2.5  High-order  weighted  compact  scheme 

In  order  to  capture  the  shock  without  elimination  of  small  vortex  which  is  particularly 
important  to  DNS  for  shock-boundary  layer  interaction,  we  developed  a  so-called  ‘high- 
order  weighted  compact  scheme’.  The  details  will  be  described  in  Chapter  3. 

2.6  Non  reflecting  boundary  conditions 

In  DNS/LES  for  compressible  flow,  especially  for  complex  geometry,  the  properly 
imposed  boundary  condition  is  required  to  prevent  the  wave  from  reflecting  at  the 
boundary  and  contaminating  the  inner  region  of  a  computational  domain.  The 
specification  of  the  outflow  boundary  conditions  becomes  one  of  the  major  difficulties. 
To  avoid  reflection  of  outgoing  waves,  a  buffer  domain  (Street  &  Macaraeg,  1989)  was 
introduced  adding  to  the  original  computational  domain  at  the  outflow  boundary.  The 
governing  equations  in  the  buffer  domain  have  been  modified  by  increasing  diffusion  in 
certain  direction.  Alternatively,  a  sponge  layer  approach  has  been  used  by  a  number  of 
authors.  Both  buffer  and  sponge  zones  are  capable  of  absorbing  the  outward  moving 
waves,  and  have  been  used  in  many  DNS/LES  computations.  But  these  methods  have 
some  shortages.  First,  extra  sponge  or  buffer  areas  have  to  be  added  to  the  original 
domain,  which  will  increase  the  number  of  grid  points  and  computational  cost.  Second, 
the  sponge  approach  can  only  be  applied  when  the  equations  of  perturbation  are 
considered,  which  may  not  be  applicable  when  the  so-called  "base  flow”  does  not  exist  or 
can  hardly  be  defined.  An  alternative  approach  to  buffer  or  sponge  method  is  the  non¬ 
reflecting  boundary  condition  based  on  the  characteristic  analysis.  The  characteristic 
boundary  conditions  for  Euler  system  can  be  found  in  a  number  of  literatures.  The 
concept  of  non-reflecting  boundary  conditions  was  proposed  by  Thompson  (1987,  1990), 
where  the  idea  of  specifying  the  boundary  conditions  according  to  the  inward  and 
outward  propagating  waves  has  been  introduced.  Usually  the  outgoing  waves  have  their 
behavior  defined  entirely  by  the  solution  at  and  within  the  boundary,  and  no  boundary 
conditions  are  needed.  Therefore,  the  number  of  boundary  conditions  is  equal  to  the 
number  of  incoming  waves.  In  Thompson's  paper,  boundary  conditions  of  Euler 
equations  are  given  for  different  situations  without  considering  the  viscosity  effect.  The 
characteristic  boundary  conditions  for  Euler  system  cannot  be  directly  used  for  Navier- 
Stokes  equations.  Actually,  there  is  no  characteristic  relation  in  Navier-Stokes  system. 
However,  in  some  circumstances,  such  as  free  stream,  far  field,  etc.  where  the  viscosity  is 
less  important,  the  characteristic  relation  can  still  be  used.  Poinsot  and  Lele  (1992)  have 
extended  Thompson's  method  to  specify  the  boundary  conditions  for  the  Navier-Stokes 
equations,  where  the  effect  of  viscous  has  been  taken  into  account.  However,  they  only 
give  the  expression  of  boundary  conditions  in  the  Cartesian  coordinates.  Based  on  the 
previous  work  by  the  above  authors,  we  have  developed  non-reflecting  boundary 
conditions  for  compressible  flow  in  curvilinear  coordinates  (Jiang  et  al,  1999b).  The 
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computation  shows  our  work  is  completely  successful  for  eliminating  non-physical  wave 
reflection.  Apparently,  it  is  an  original  work  and  important  to  DNS  for  curvilinear 
coordinates.  Based  on  the  characteristic  analysis,  the  hyperbolic  terms  in  £  direction  can 
be  modified  as: 


9/?  3/?  3 u  3v  3w.  dp  3m  -  3v  „  3 w. 

_+,l+v_+^_+,_+,_)^_+^_+f  _+C_)+ra,=0 

at  drj  p  dij  p  d£ 

3v  3v  1  dp  3v  1  „  3 (p  , 

— +  d3+y— +— n  -£-  +  W— +  -r-^  +  m3  =0  (2.11) 

3/  drj  p  y  dr]  dC  P  dC 

3w  T.3w  1  dp  ...dw  l  y  dp  .  . 

37  +  ^4+^3-  +— ^3-  +  ^  — ^  — +  vis4  =0 
dt  drj  p  drj  3  £  p  d£ 

dp  dp  du  dv  dw.  ...dp  du  „  dv  „  3w. 

_+rfs+v__+^_+^_+,_,+H,_+^_+c_+c_)+ros=o 


where  visl  -  vis5  represent  viscous  terms  in  curvilinear  coordinates,  and 


*x 

d2 

d, 

d, 

d5 


l-d^+L^  +  L,] 

c  2 

WA) 

|(A+A) 


(2.12) 


In  (2.12)  c  is  the  sound  wave  speed  and  (5  =  +  <f,2  +  ^,2  . 

^represent  the  amplitude  variations  of  the  characteristic  wanes  corresponding  to  the 

characteristic  velocities,  which  are  given  by 

Al=U-C( 

A2=Ai=A4=U 
As  =U  +  CS 
where  =  c/3  and 
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2  dp  dps 


2  3<f  3<T 

T  rr,  e  9m  ~  3w 

^=^("^3?  +  ^3#)  (2-13) 

r  ....  e  3m  •„  9vw 

L‘=u^rp^ 

T  ,TT  sri oc  ,e  3m  e  3v  e  9w.  3p. 

5  ~~  eX-0&  3^  +  ^  9£  +  ^  3#  +  3£ 

These  equations  will  be  used  for  neighbors  of  boundary  points  in  £  direction.  The 
equations  for  7]  and  C,  directions  are  similar.  In  this  way,  the  non-physical  wave 
reflection  can  be  effectively  eliminated. 

2.7  Effective  sub-grid  models  (Structure-function  subgrid  model) 

Since  LES  only  resolved  the  large  length  scales  but  filtered  small  length  scales,  to  adopt  a 
right  sub-grid  model  to  represent  these  filtered  small  length  scales  becomes  particularly 
important. 

The  filtered  stmcture-function  model  is  developed  by  Ducros  et  al.(1996).  The  subgrid 
scale  shear  stress  and  heat  flux  can  be  modeled  as 


fe-JfO  2du 


dx,  dxk  J  3  dx„ 


(2.14) 


„  _  ypv,  3 f 

Qk  ~  D  -j 

Pr,  dxk 

Here  Prtis  the  turbulent  Prandtl  number  taken  equal  to  0.6  as  in  isotropic  turbulence,  v(is 
the  turbulent  kinetic  viscosity  defined  as 


v,(x,t)  =  0.0014C‘3/2A[F2(3)(x,0] 


(2.15) 


where  F23)  is  the  filtered  structure  function.  In  the  case  of  flat-plate  boundary  layer  flows 

with  meshes  flattened  in  the  wall-normal  direction,  F23)  takes  the  four-neighbor 
formulation  proposed  by  Normand  and  Lesieur  (1992). 

F(3)  =-["|!m(3)  —  ;7(3)l|2  +ll;7(3)  -;7(3)I!2  +!I/7(3)  -  ,7<3)l|24- Ilr7<3>  _,7(3)||2" 

2  4|_lr«+».y  uu  1  +Ir«-W  uu  I  +  ||Mi.;+i  ui,j  I  +  \\ui,j-i  ui.j  I 


where 


Ki  =HPw(uu) 
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A  =  (AfAy)1/2  is  used  to  characterize  the  grid  size,  where  Axand  Ay  denote  the  grid  size 

in  the  streamwise  direction  and  spanwise  direction,  respectively,  C*  is  the  Kolmogorov 
constant  taking  the  value  of  1.4.  HP(3)  is  a  discrete  Laplacian  filter  iterated  3  times,  which 
is  served  as  a  high-pass  filter  before  computing  the  structure  function.  The  first  iteration 
of  the  Laplacian  filter  HF*l)  is  defined  by 


u«l=HPw(uiJ)  =  ui+lJ-2uiJ+u, 


1,7  +  «i  J+I 


■2“ ij+Z'j-i 


2.8  MPI  parallel  computation 


(2.16) 


Message  Passing  Interface  (MPI,  e.g,,  Pacheco,  1994),  which  is  a  standard  message¬ 
passing  library  for  writing  message-passing  style  parallel  programs,  provides  us  a  better 
route  to  write  parallel  programs.  MPI  consists  of  a  small  group  of  functions  that  are  used 
to  support  communication  among  processors.  A  message-passing  function  is  simply  a 
function  that  explicitly  controls  data  transmitting  from  one  processor  to  the  other. 
Because  MPI  has  become  a  standard  library,  most  existing  parallel  computers  support  it. 
Therefore,  an  MPI  programmer  has  no  need  to  worry  about  portability.  These  computers 
for  MPI  include  not  only  distributed-memory  parallel  computers  but  also  shared-memory 
parallel  computers.  Even  on  workstations,  MPI  can  also  be  used.  This  greatly  reduces  the 
limitation  of  parallel  application  in  DNS. 

In  our  previous  effort,  an  MPI  parallel  code  was  developed  based  on  our  serial  code.  The 
parallelism  of  the  DNS  code  is  examined.  The  performance  of  the  parallel  program  is 
examined  for  our  compressible  DNS  code  on  an  SGI  Origin  2000  computer  and  the 
results  show  great  computing  efficiency  of  parallel  machines  with  MPI  (see  Figure  3) 


Figure  3.  Linearly  Speed  up  of  MPI  Computation 


In  order  to  maintain  good  balance  data  must  be  uniformly  distributed  among  multiple 
processors.  Here,  because  the  computation  is  related  to  each  grid  point,  how  to  distribute 
these  points  among  processors,  or  how  to  partition  the  computational  domain  is  very 
important  to  improve  parallelism.  Since  MPI  code  is  very  portable,  our  LES  code  can  be 
used  for  any  parallel  machine  including  our  own  machine  and  DOD  supercomputers. 


16 


Chapter  3 

High-Order  Weighted  Compact  Scheme 

3.1  High-order  Compact  Scheme 

Recently  compact  schemes  have  been  widely  used  in  CFD.  Standard  finite  difference 
schemes  need  to  be  at  least  one  point  wider  than  the  desired  approximation  order.  It  is 
difficult  to  find  suitable  and  stable  boundary  closure  for  high  order  schemes.  However, 
the  compact  scheme  can  achieve  higher  order  without  increasing  the  stencil  width.  As 
the  compact  scheme  has  an  implicit  form  and  involves  neighboring  grid  point  derivative 
values,  additional  free  parameters  can  be  used  not  only  to  improve  the  accuracy  but  also 
to  optimize  the  other  properties  such  as  resolution,  stability,  and  conservation.  A  family 
of  centered  compact  schemes  proposed  by  LeLe(Lele,  1992)  have  been  proved  to  have 
spectral-like  resolution.  The  conservation  property  is  also  important,  especially  for  flow 
with  shocks. 

Though  the  advantages  of  compact  schemes  are  obvious,  there  are  still  difficulties  in 
using  them  to  solve  problems  involving  shock  waves  or  discontinuities.  When  they  are 
used  to  differentiate  a  discontinuous  function,  the  computed  derivative  has  grid  to  grid 
oscillation.  Recently  the  ENO  (Harten,  1987;  Shu,  1988,  1989)  and  WENO  (Liu,  1994; 
Jiang,  1996)  schemes  have  been  widely  used  for  shock  wave  capturing  and  have  been 
proved  very  successful.  These  schemes  check  the  smoothness  of  the  candidate  stencils. 
The  ENO  scheme  selects  the  smoothest  stencil,  while  the  WENO  scheme  uses  all  the 
candidate  stencils  but  with  assigned  weights.  Inspired  by  the  success  of  the  WENO 
scheme,  we  have  developed  a  new  compact  scheme  (Jiang  et  al,  2001a)  so  that  the  new 
compact  scheme  not  only  preserves  the  properties  of  compact  schemes  but  also  can  be 
used  for  shock  wave  capturing.  This  new  scheme  preserves  the  characteristic  of  standard 
compact  schemes  achieving  high  order  accuracy  and  high  resolution  by  a  compact  stencil. 
The  improvement  of  this  new  scheme  over  the  standard  compact  scheme  is  that  it  can 
accurately  capture  shock  waves  without  oscillation.  The  idea  of  the  Weighted  Compact 
Scheme  is  similar  to  the  WENO  scheme  (Jiang  and  Shu  1996).  In  the  WENO  scheme, 
each  of  the  candidate  stencils  is  assigned  a  weight  that  determines  the  contribution  of  this 
stencil  to  the  final  approximation  of  the  numerical  flux.  The  weights  are  defined  in  such 
a  way  that  in  smooth  regions  it  approaches  certain  optimal  weights  to  achieve  a  higher 
order  of  accuracy,  while  in  regions  close  to  discontinuity,  the  stencils  that  contain  the 
discontinuities  are  assigned  a  nearly  zero  weight.  Similarly  to  this  idea,  the  Weighted 
Compact  Scheme  is  constructed  by  the  combination  of  the  approximations  of  derivatives 
on  candidate  stencils.  The  stencil  that  contains  the  discontinuity  has  less  contribution.  In 
this  way,  the  oscillation  near  the  discontinuity  can  be  reduced,  while  the  characteristics  of 
the  compact  scheme  can  still  be  preserved.  The  building  blocks  of  the  Weighted  Compact 
Scheme  are  the  standard  compact  schemes.  The  Weighted  Compact  Scheme  is  a  hybrid 
of  different  forms  of  standard  schemes.  Compared  to  the  WENO  scheme,  this  Weighted 
Compact  Scheme  can  achieve  higher  order  accuracy  and  higher  resolution  with  the  same 
stencil. 
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3.2  High-order  Weighted  compact  scheme 


For  simplicity,  we  consider  a  uniformly  spaced  grid  where  the  nodes  are  indexed  by  j. 
The  independent  variable  at  the  nodes  is  xj  =  h(j  - 1)  for  1  <  j  <N  and  the  function 

values  at  the  nodes  fJ  =  f[xj)  are  given.  The  finite  difference  approximation  /'  to  the 


first  derivative  of  the  function  /  on  the  nodes  can  be  written  in  the  following  general 
form  while  the  standard  compact  scheme  (Lele,  1992)  is  used. 


P-/^  +aj;_l +/;+«,/;,  +p+/;2  =-(bjj_2+a_fM  +c/.  +ajJ+i  +bjm) 


For  the  point  j,  we  define  three  candidate  stencils  containing  the  point  j : 


(3.1) 


So  [xj-2  >  xj-i  >  xj )’  S,  [xj~\  >  Xj .  Xj+i )>  S2  [xj ,  xj+]  ,xj+2 ) 


On  each  of  them  we  can  get  a  compact  scheme.  By  matching  the  Taylor  series 
coefficients  to  various  orders,  the  third  and  fourth  order  compact  schemes  corresponding 
to  each  stencil  can  be  derived.  The  coefficients  are  given  as  follows: 


S0:  P_  =  P,  a_=2p  +  2,  b_  a_  =  2fi-2, 

c  .  _  1  1  3  3 

51  .  a_  ,  a+  —  ,  a_  —  ,  a+  —  ,  c  —  0, 

4  4  4  4 

52  :  P+  =  P,  a+=2P  +  2,  b+  =|p+|,  a+=-2P  +  2, 


c  = 


5 

2’ 


where  P  is  a  free  parameter.  The  coefficients  that  are  not  listed  are  set  to  zero.  The 
schemes  corresponding  to  stencils  S0  and  S2  are  third  order  one-sided  schemes,  and  the 

scheme  corresponding  to  Sx  is  a  fourth  order  centered  scheme.  With  these  three  sets  of 
coefficients,  we  get  three  different  equations  from  Eq.  (3.1).  These  equations  are 
represented  by  F&  Fj,  F2.  When  these  equations  are  assigned  specific  weights  and 
combined,  a  new  scheme  is  obtained: 


f  =  c0f0  +  c,f1+c2f2 

where,  C0  +  Q  +  C2  =  1 .  If  the  weights  are  properly  chosen,  the  new  scheme  can  achieve 
a  higher  order  because  the  additional  free  parameters  are  introduced.  If  we  set 


•  =  r  -  1  r  -8~12P 

0  2  18-24p’  1  9-12P 


(3.2) 


the  new  scheme  is  at  least  a  sixth  order  centered  compact  scheme.  This  process  implies 
that  this  sixth  order  centered  compact  scheme  can  be  represented  by  a  combination  of 
three  lower  order  schemes.  Obviously,  the  scheme  F  is  a  standard  compact  scheme  and 
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cannot  avoid  the  oscillation  near  the  discontinuity.  Can  we  define  the  weights  in  such  a 
way  that  the  scheme  has  the  non-oscillatory  property?  Then  the  idea  of  the  WENO 
scheme  is  introduced  to  determine  the  new  weight  for  each  stencil.  The  weights  are 
determined  according  to  the  smoothness  of  the  function  on  each  stencil.  According  to  the 
WENO  method,  the  new  weight  (ak  is  defined  as  (Jiang  and  Shu,  1996): 


Jt=0 


(3.3) 


where,  e  is  a  positive  small  number  that  is  used  to  avoid  the  denominator  becoming  zero. 
The  smoothness  measurement  ISk  is  defined  as  following: 


K.  =  2fj-,  +  /,)*  +±(/>,  -4  +3  f,f 

IS,  =  || ~2fj  +/,,,)=  +  ~f,J  (3.4) 

</,  -2/;.,  +  /,„)’  “ 4/;>,  +3/;)2 

In  smooth  regions,  the  following  forms  can  be  obtained  by  a  Taylor  expansion  of  (4) 

/5o=ll(/V)  +  i(2//*"l/^3)  +0(^6) 

7S,=l|(/^2)  +  ^(2/^  +  I/<r*3)  +°(ft6)  (3'5) 

^2=^(/V)+|f2rA-|/^0  +o(a6) 

12  4\  3  J 

If  /'  =£  0 ,  then 

is0=(,f'hf  l  +  +0(h6) 

lSx={f'h)2  1+  /f- -  +  V  ^  +0(fc6) 

\ 12/  o  J 


Through  a  Taylor  expansion,  it  can  be  easily  proved  that  the  new  weight  03*  satisfies 
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Obviously,  oo0  -  co2  =  0(h3) ;  this  ensures  that  the  new  scheme 


F  =  <o0F0  +  tOjFj  +  g>2F2 


(3.6) 


still  can  be  of  sixth  order  .If  /'  =  0 ,  it  can  also  be  proved  that  the  new  scheme  (3.6)  is 
of  sixth  order.  The  conclusion  is  that  the  Weighted  Compact  Scheme  (3.6)  can  preserve 
the  higher  order  accuracy  in  the  smooth  region  when  the  WENO  weights  are  used. 
According  to  these  definitions,  in  the  smooth  region  the  new  scheme  is  nearly  a  standard 
sixth  order  centered  compact  scheme.  The  non-smooth  stencil  is  assigned  a  small  weight 
so  that  the  non-oscillatory  property  is  achieved.  The  coefficients  of  the  final  Weighted 
Compact  Scheme  are  given  as  follows: 


P-  =N0,  a_  =(2p  +  2)co0+ito1,  a+ =(2P  +  2)o)2 +^0),,  p+=P(o2, 


-VI 

2  2 


,  o_=(2|3-2K-|o>1.  «=ge+f)«b-gp+£ 

a,=(-2P  +  2K+!»„  4*=(fP  +  l)“2 


(0 


2 5 


(3.7) 


If  P  =  0,  the  scheme  is  tridiagonal;  otherwise,  it  is  pentadi agonal.  The  free  parameter  P 
can  be  used  to  optimize  the  scheme  when  the  properties  of  high  resolution,  conservation, 
and  stability  are  concerned.  From  the  above,  we  can  find  that  the  weights  play  a  very 
important  role  in  the  Weighted  Compact  Scheme.  They  can  be  used  to  optimize  the 
accuracy.  In  addition,  these  weights  can  make  the  new  scheme  non-oscillatory.  In  this 
proposed  project. 

3.3  Computational  Results  for  rectangular  grids 
1)  Convection  Equation 

We  first  solve  the  following  model  equation  with  different  initial  functions 
«,  +1^=0  -1  <  ;c  <  1 

u{x, 0)  =  w0  (x)  periodic  with  a  period  of  2  (3-8) 

The  initial  functions  are 
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N 

Li  error 

Li  order 

N 

Li  error 

Li  order 

20 

2.03D-4 

20 

1.29D-2 

40 

3.01D-6 

6.07 

40 

2.11D-4 

5.93 

80 

4.59D-8 

6.04 

80 

3.15D-6 

6.06 

160 

6.17D-10 

6.23 

160 

4.32D-8 

6.19 

Table  1 

Table  2 

Tables  1  and  2  list  the  errors  and  orders  of  accuracy  in  the  Lj  norm  between  the 
numerical  result  and  the  exact  solution  at  time  t  =  1  for  initial  functions  (1)  and  (2), 
respectively.  N  is  the  number  of  grid  points.  High  order  accuracy  for  these  two  smooth 
initial  functions  has  been  achieved. 

Figures4-6  illustrate  the  results  at  t  =  0.5  for  initial  functions  (3),  (4),  and  (5).  Clearly 
there  is  no  numerical  oscillation  observed  in  the  region  near  the  discontinuity. 


Figure  4.  The  solution  at  t  =  0.5  FigureS.The  solution  at  t  =  0.5  Figure  6.  The  solution  at  t  =  0.5 
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2)  Burgers’  Equation 


The  Weighted  Compact  Scheme  is  applied  to  the  nonlinear  Burgers’  equation.  With  the 
given  initial  condition,  the  exact  solution  will  develop  a  moving  shock  wave. 


-!<*<! 


u(x,0)  =  —  +  sin(7ix)  periodic  with  a  period  of  2 


Fig.7  (a)  and  (b)  show  the  wave  at  different  times.  The  shock  appears  at  t  =  0.55  and  is 
accurately  captured  by  the  scheme.  No  obvious  oscillation  is  observed. 


Figure  7.  The  solution  of  Burgers’  equation  at  (a )t  =  0.3183  and  (b)  f  =  0.55 


3)  Euler  Equation 

We  have  applied  the  scheme  to  the  Euler  equation.  The  equations  are: 

dU  dF  . 

T~  +  T"  =  0 
at  ox 

U  -  (p,  pu,  E)t  ,  F  =  (pu,pu  +p,u(E+p)J 
The  initial  conditions  are: 

v  JM.i)  *<° 

0  [(0.125,0,0.1)  x>0 

The  results  are  shown  in  Fig.  8. 


(3.9) 
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(a)  (b) 

Figure  8.  The  shock  wave  solutions  of  Euler  equation  at  t  =  2 .  Number  of  points  =  100, 
From  these  preliminary  applications  of  the  new  Weighted  Compact  Scheme,  we  can  see 
the  prospects  of  this  new  scheme.  When  the  scheme  is  optimized  for  optimal  resolution 
and  conservation,  the  shock-capturing  property  should  be  much  better. 


3.4  Computational  results  for  2-d  curvilinear  grids 

For  2-d  curvilinear  coordinates,  several  testing  cases  have  been  conducted.  The  results 
show  that  the  Weighted  Compacted  Scheme  we  developed  (Jiang  et  al,  2001a)  can  be 
used  for  curvilinear  coordinate  system.  The  medium  distortion  of  the  grid  doesn’t  have 
much  effect  on  the  solutions. 

1)  Shock-fluctuation  interaction  problem 

2D  Euler  equations  are  solved  with  the  following  initial  conditions: 


23 


pre  -  shock  :x>l  w,  =  -c,  sin#  cos(;d  cos  0+yk  sin  6) 
v,  =  — c,  cos9cos(xk  cos  0  +  yk  sin  9) 


Pi  =1 


post  -  shock  :  x  <  1  h. 


2(M 2  - 1) 
(/+1)M2' 


v2  =  0 

(y  +  l)M2 

p2=(l+^ — ’-)PX 

y  + 1 

where  M=8,  which  indicates  a  very  strong  shock  wave.  The  following  Figure  9  show  the 
results  obtained  on  different  grids.  It  can  be  seen  that  the  sharpness  of  the  shock  is  kept 
very  well  even  the  grid  is  distorted. 


Grid=60X40 


2)  Vortex  pairing  in  a  mixing  layer 

Vortex  pairing  process  is  forced  by  adding  velocity  disturbances  to  the  initial  mean 
velocity  profile  with  opposite  free-stream. 

Mean  flow  and  temperature  fields  are  given  by 
u  =  tanh(2y) 

f  =  l  +  -^M2(l-«2) 

Disturbances  are  added  as 
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where  Ax  =  0.05,  A2  =  0.025,  B  =  10,  Lx  =  20,  Ly  =  40,  Re  =  2000,  M  =  0.8 . 

The  following  Figure  10  show  the  appearance  of  the  shock  waves  and  their 
interactions  with  the  vortices.  The  results  show  that  the  position  and  the  shapes  of  the 
shock  waves  are  not  affect  by  the  distortion  of  the  grid. 


t-  +  +  Grid  +  +  +  +  +  +Votfic'rty  +  + 


+  +  ^Pressure  + 


Figure  10.  Comparison  between  rectangular  and  curvilinear  grids  for  vortex  pairing 


•+ 


Conclusion 


The  weighted  compact  scheme  can  preserve  the  high  order  accuracy  in  the  smooth  area 
while  has  no  oscillation  for  discontinuity  (shocks).  The  scheme  works  well  for  both 
rectangular  and  curvilinear  grids. 
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Chapter  4 


DNS  for  Flow  Separation  and  Transition  around  Airfoils  with  Attack 
Angles 

Flow  transition  in  separation  bubbles  is  a  classic  topic  and  has  been  studied  for  many 
years  (Boiko  et  al,  2002).  However,  most  of  work  was  focused  on  flow  around  a  hump 
placed  on  a  flat  plate  (Musad  et  al,  1994)  or  for  a  blunt  leading  edge  (Yang  &  Voke, 
2001).  Flow  separation  and  transition  around  an  airfoil  with  attacked  angle  is  rarely  found 
in  literature  due  to  its  complexity.  The  linear  stability  theory  (LST;  see  Drazin  &  Reid, 
1981)  is  mainly  a  local  analysis  with  assumption  of  parallel  base  flow.  The  parabolized 
stability  equations  (PSE;  see  Bertolotti,  1992)  assume  a  steady  base  flow  with  no  elliptic 
part.  These  assumptions  do  not  apply  for  the  case  of  flow  separation  and  transition 
around  airfoil  with  attack  angle  where  no  steady  base  flow  exists  and  the  transition 
process  is  dominated  by  an  elliptic  process  especially  in  the  late  stages.  Though  it  is  true 
that  the  LST  and  PSE  cannot  provide  a  correct  prediction  for  the  case,  the  Kelvin- 
Helmoholtz  instability  mechanism  in  a  separated  shear  layer,  which  is  obtained  by 
inviscid  stability  theory,  still  dominates  in  the  early  transition  stage. 

Boundary  layer  separation  and  transition  exist  in  many  engineering  flows  around 
wings  and  blades.  When  an  adverse  pressure  on  a  laminar  boundary  layer  over  a  surface 
is  strong  enough,  the  laminar  boundary  layer  separates  from  the  surface.  Separation  and 
transition  in  these  types  of  flows  are  strongly  coupled.  The  instability  at  the  separation 
zone  is  widely  accepted  as  dedicated  by  the  Kelvin-Helmholtz  mechanism.  Transition 
takes  place  owing  to  nonlinear  breakdown  of  spatially  growing  traveling  waves  in  the 
separated  free  shear  layer  (Yang  &  Voke,  2001).  When  the  shear  layer  becomes  turbulent, 
the  detached  shear  layer  may  reattach  to  the  surface,  creating  a  separation  bubble  and 
forming  attached  turbulent  boundary  layer.  Obviously,  the  length  of  the  bubble  is  closely 
related  to  when  and  where  the  transition  takes  place.  On  the  other  hand,  the  size  of  the 
bubble  could  directly  affect  the  flight  characteristics  of  the  airfoil  and  the  efficiency  of 
the  turbine  machine.  The  flow  separation  over  a  wing  in  flight  will  cause  loss  of  the  lift 
and  increase  of  the  drag  which  threatens  the  stability  and  efficiency  of  the  aircraft. 
Understanding  the  mechanism  of  the  separation  and  transition  is  of  great  importance  in 
improving  design  in  aircraft  and  turbo-machinery. 

With  the  development  of  computer  resources  and  efficient  numerical  methods,  high 
resolution  and  high  accuracy  direct  numerical  simulation  has  been  becoming  feasible  on 
flow  separation  and  transition  around  airfoil  with  low  Reynolds  number.  In  this  work, 
high-accuracy  and  high-resolution  numerical  simulations  are  conducted  to  investigate  the 
details  of  separation  and  transition  processes  on  a  NACA  0012  airfoil  The  organization  of 
the  paper  is  described  as  follows.  In  section  1  and  2,  the  governing  equations  and 
numerical  methods  will  be  introduced  briefly.  In  section  3,  the  definition  of  the  study 
case  and  the  computational  set  up  will  be  given.  In  section  4,  the  computational  results 
will  be  presented  and  analyzed. 
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4.1  Governing  Equations 


The  three-dimensional  compressible  Navier-Stokes  equations  in  generalized 
curvilinear  coordinates  (<f,  r),  ip)  are  written  in  conservative  forms: 

IM  +  M  +  M  +  M  =  0  (4.1, 

J  dt  dg  dp  d£ 

The  vector  of  conserved  quantities  Q,  inviscid  flux  vector  (E,  F,  G),  and  viscous  flux 
vector  ( Ev ,  FV,GV)  are  defined  via 


pu 

Q  =  pv 


pU 

pUu  +  pgx 
pUv+p^y 
pUw+p £ 
Me,  +  p) 


pv 

pVu  +  prjx 


pVw+  prjz 
J(E,  +  p) 


pw 

pWu  +  pCx 


pVv  +  PVy  ,G=-  pWv  +  pC 


pWw  +  p£z 

ME,  +  p) 


t  E  +r  £  +r  t 

xx^x  yx^y  zx^>  z 


E*=~ 

+  *Jt 

kQX^X  +  Qy^y  +  QZ£,Z 


t  71  +  T  7]  +  T  71 

xx  ix  yx  / y  zx  tz 


K=j  rxyVx+ryyrjy  +  rzypz 
VxzVx+VyzVy+'CJh 

AVx+Qvr?Y+Qzrjz 


TxxCx  +  TyxCy  +  TuCz 
^V~~J  ^yy^y  ^ zy^ z 

^xzCx  ^  yz C y  +  ^zzCz 

^QxCx+QyCy  +QZCzy 

'x/  p  y\ 

where  J  =  ■  .  — — r-  is  Jacobian  of  the  coordinate  transformation  between  the  curvilinear 
o{x,y,z) 

(£  V>  C)  and  Cartesian  {x,y,z)  frames,  and  ^x,qy,^z,px,py,pz,CxXy^z  are  coordinate 
transformation  metrics.  The  contravariant  velocity  components  U,V,W  are  defined  as 
U=ugx+vgy  +w£,  V  =  urjx  +vrjy  +wl]z ,  W  =  uCx+vCy+wCz. 
Qx’Qy’Qz  in  energy  equation  are  defined  as 

Qx  =-^lx+UTxx+VTxy+W^xz 
Qy  =  ~q,  +  ur^  +  VT^  +  wryz 
Qz  =~qz  +uTxz+vTyz+wrzz 

Et  is  the  total  energy.  The  components  of  the  viscous  stress  tensor  and  heat  flux  are 
denoted  by  rxx,tyy,ta,Txy  ,rxz,ryz,  and  qx,qy,qz,  respectively. 


In  the  dimensionless  form,  the  reference  values  for  length,  density,  velocities, 
temperature,  pressure  and  time  are  L,pr,Ur,Tr,  prU2 ,  and  LIU r  respectively.  In  this 
study,  the  free  stream  parameters  are  chosen  as  reference  values.  The  chord  length  of  the 
airfoil  is  used  as  the  reference  length.  The  dimensionless  parameters  arise  from  non- 
dimensional  Mach  number  M r ,  Reynolds  number  Re ,  Prandtl  number  Pr ,  and  the  ratio  of 
specific  heats  y ,  are  defined  respectively  as  follows 


~M-  = 


U, 


4W 


Re  = 


Pr 


Pr  = 


CpPr 


y  =  . 


where  R  is  the  ideal  gas  constant,  Cp  and  Cv  are  specific  heats  at  constant  pressure  and 
constant  volume,  respectively.  Through  out  this  work,  Pr  =  0.7  andy  =  1.4 .  Viscosity  is 
determined  according  to  the  Sutherland’s  law  in  dimensionless  form 


P  = 


tV2{i+s)  c_ 

T  +  S 


110.3K 

r„ 


The  governing  system  is  closed  by  the  equation  of  state. 
jM2p  =  pT 

E  =  —  +  —  p{u2  +  v2  +w2) 
y- 1  2  x  } 

The  components  of  the  viscous  stress  tensor  and  heat  flux  in  non-dimensional  form  are  as 
follows: 


T..=JL 

11  Re 


A 


Qi  = 


dui  du 
KdXj  dxf  j 

M 


3  y 


dT 


(y-l)M2RePr  dxl 


dxk 


4.2  Numerical  Methods 


The  numerical  method  used  in  this  study  has  high  order  accuracy  and  high 
resolution.  The  governing  equations  are  solved  using  LU-SGS  implicit  method  based  on 
second  order  Euler  Backward  scheme.  The  sixth-order  centered  compact  difference 
scheme  is  used  for  spatial  derivatives.  The  eighth-order  compact  filter  is  employed  to 
reduce  numerical  oscillation. 

In  Eq.  (1),  a  second  order  Euler  Backward  scheme  is  used  for  time  derivatives,  and 
the  fully  implicit  form  of  the  discretized  equations  is  given  by 
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Qn+l  is  estimated  iteratively  as: 
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where. 


Qn+1=QP+SQP 
SQP  =  Qp+l  -Qp 
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At  step  p  =  0 ,  Qp  =  Qn ;  as  SQP  is  driven  to  zero,  Qp  approaches  Qn+t .  The  flux  vectors 
are  linearized  as  follows: 


En+l  ~EP+APSQP 
F"+1  «  Fp  +BPSQP 
Gn+1  ~GP  +CPSQP 


So  that  Eq.  (4.2)  can  be  written  as: 
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SQP  =  R 


(4.3) 


where  R  is  the  residual: 

*  =  -(|e' -2Q-  +iQ’"']-4(j[Di(E-£,)+D,(F-F.)  +  Df(G-Gv)]'(4.4) 

Dg,Dv,D(  represent  partial  differential  operators,  and  A,  B,  C  are  the  Jacobian  matrices 
of  flux  vectors: 

A-®  bJP  g  =  ~ 
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The  right  hand  side  of  Eq,  (4.3)  is  discretized  using  sixth-order  compact  scheme  (Lele, 
1992)  for  spatial  derivatives,  and  the  left  hand  side  of  the  equation  is  discretized 
following  LU-SGS  method  (Yoon,  et  al .,  1992)  In  this  method,  the  Jacobian  matrices  of 
flux  vectors  are  split  as: 

A  -  A+  +  A~ ,  B  =  B++B~,  C  =  C++C" 

where 


A*=^U±r,/].  C*=|[C±rc/] 

and, 

rA  =  armax|>l(A)|]+  v ,  rB=  jnnax(^(5)|]+  v ,  rc  =  jcmax j^(c]|]+  v 

where  /1(a),  A{b),  A(c)  are  eigenvalues  of  A,  B,  C  respectively,  /cis  a  constant  greater 

than  1.  v  is  taken  into  account  for  the  effects  of  viscous  terms,  and  the  following 
expression  is  used: 
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The  first-order  upwind  finite  difference  scheme  is  used  for  the  split  flux  terms  on  the  left 
hand  side  of  Eq.  (4.3).  This  does  not  affect  the  accuracy  of  the  scheme  when  solutions  are 
converged.  As  the  left  hand  side  is  driven  to  zero,  the  discretization  error  will  also  be 
driven  to  zero.  The  finite  difference  representation  of  Eq.  (4.3)  can  be  written  as: 
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In  LU-SGS  method,  Eq.  (4.5)  is  solved  with  three  steps.  First  initialize  SQ°  using 
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In  the  second  step,  the  following  relation  is  used: 
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For  the  last  step,  $2P  is  obtained  by 


■I  +  Atj(rA+rt+rc)l 


-1 

Atj[A‘^2, 


p  +7? 


■«S2,W+C" 


*+l 


The  sweeping  of  the  computational  domain  is  performed  along  the  planes  of 
i  +  j  +  k  =  const ,  i.e.  in  the  second  step,  sweeping  is  from  the  low-left  comer  of  the  grid 
to  the  upper-right  comer,  and  then  vice  versa  in  the  third. 

Parallel  computation  based  on  Message  Passing  Interface  has  been  utilized  to 
improve  the  performance  of  the  code.  The  parallel  computation  is  combined  with  the 
domain  decomposition  method.  The  computational  domain  is  divided  into  n  equal-sized 
sub-domains  along  £,  direction  which  is  usually  the  streamwise  direction.  The  details  of 
the  numerical  methods  and  parallel  computation  can  be  referred  to  our  previous  work 
(Jiang,  et  al,  1999a;  Shan,  et  al,  2001). 


4.3  Problem  Definitions  and  Boundary  Conditions 


Numerical  simulations  are  performed  on  a  NACA0012  airfoil  at  attack  angle  of  4°. 
Free  stream  velocity  U„,  pressure  pm,  temperature  T„  and  chord  length  of  the  airfoil  C 
are  selected  as  reference  velocity,  pressure,  temperature  and  length  respectively,  which 
are  used  to  non-dimensional  governing  equations.  The  computational  domain  is  plotted  in 
Figure  11.  The  upstream  boundary  is  3  chord  lengths  away  from  the  leading  edge  of  the 
airfoil.  The  upper  and  lower  boundaries  are  about  4  chord  length  from  the  solid  surface. 
The  outflow  boundary  is  2  chord  lengths  downstream  of  the  trailing  edge.  The  airfoil  is 
regarded  as  infinitely  long  in  the  spanwise  direction.  In  our  simulations,  spanwise  length 
is  set  as  Ly  =  0.1C  and  periodic  boundary  condition  is  imposed  at  the  spanwise 
boundaries. 
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Figure  11.  Computational  domain 

The  flow  and  computational  conditions  are  listed  in  Table  3.  The  Reynolds 
Number  based  on  free  stream  velocity  and  chord  length  isRe  =  105 .  Free  stream  Mach 
Number  is  M „  =0.2.  The  numbers  of  grid  point  in  £ ,  t] ,  £  directions  are  Nx  =1200, 
Ny  =  32 ,  Nz  =180  respectively.  Grid  distributions  in  the  (x,  z)  plane  and  on  the  airfoil 

surface  are  shown  in  Figure  12.  Grids  are  uniform  in  the  spanwise  direction  and  stretched 
in  the  wall-normal  direction.  C-grid  is  used  in  the  (x,  z)  plane.  The  grid  sizes  in  wall  unit 
are  also  listed  in  Table  3.  Parallel  computing  is  based  on  domain  decomposition.  The 
computational  domain  is  divided  evenly  into  N  sub-domains  along  £  direction  when  N 
processors  are  used.  In  this  work,  24  processors  are  used. 


Re=t/_C/y„ 

AOA 

NxNxNr 

*  y  z 

Ax+ 

Ay+ 

Az+ 

105 

0.2 

4° 

1200x32x180 

<13 

<15 

<1 

Table  3.  Flow  and  computational  conditions 

¥ 
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(a)  Grid  in  (x,  z)  plane 


(b)  Grid  on  the  airfoil  surface 


Figure  12.  Grid  distribution  (one  out  of  three  grid  point  is  shown) 
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As  the  flow  is  subsonic  in  free  stream,  the  uniform  free  stream  velocities  and 
temperature  are  prescribed  at  upstream  and  far  field  boundaries.  Then  density  is  decided 
by  non-reflecting  boundary  condition  developed  in  our  previous  work  (Jiang,  et  al, 
1999a).  Non-slip,  adiabatic  boundary  condition  is  used  on  the  surface  of  the  airfoil.  Non¬ 
reflecting  boundary  condition  is  applied  to  the  outflow  boundary. 

4.4  Results  and  Analysis 

In  this  study,  simulation  is  carried  out  with  a  time  step  equal  to8.35xlO“5C/t/„, . 
The  corresponding  CFL  number  is  around  400.  2D  solutions  are  used  as  the  initial  field. 
The  time  integration  has  reached  t  =  3.474C/[/„  .  Mean  values  are  obtained  by 
performing  averaging  in  the  spanwise  direction  and  in  time  over  a  period  of3C/£/„ . 

4.4.1  2D  instability  waves  and  vortex  shedding 

2D  instability  waves  and  vortex  shedding  appear  in  both  2D  and  3D  simulations. 
The  generation  and  growing  of  2D  instability  waves  is  the  initial  stage  of  the  transition. 
In  2D  simulation,  because  of  the  lack  of  3D  motions,  we  can  not  see  the  late  stage  of 
transition  and  breakdown.  But  the  2D  simulation  can  provide  us  important  information  on 
how  the  instability  waves  develop  and  grow.  2D  simulation  starts  from  a  uniform  flow 
field.  This  initial  field  is  not  the  solution  of  the  governing  equations  and  may  bring  in 
some  disturbances  (initial  disturbance  brought  by  numerical  initial  condition).  When 
perturbations  introduced  by  the  initial  field  can  not  dissipate  through  the  time  integration, 
they  may  trigger  the  most  unstable  instability  waves.  This  is  the  reason  why  we  see 
unsteady  behaviors  in  the  simulation  although  all  boundary  conditions  are  steady  and  no 
external  disturbances  are  enforced.  Flow  separation  and  vortex  shedding  appear  on  the 
suction  surface  of  the  airfoil  (Figure  13),  where  a  separated  mixing  layer  and  vortex 
shedding  are  clearly  demonstrated  by  plotting  contours  of  instantaneous  spanwise 
vorticity.  The  separation  zone  can  be  seen  clearly  from  the  time  averaged  velocity  vectors 
shown  in  Figure  14,  but  no  re-attachment  is  observed  for  2-D  simulation  (2-D  simulation 
cannot  represent  the  real  physics).  The  2D  solutions  are  well  developed  and  quasi- 
periodic  behavior  is  built  up.  It  can  be  seen  clearly  that  the  boundary  layer  separates  from 
the  airfoil  surface  and  develops  to  unstable  shear  layer  which  leads  to  vortex  shedding. 
These  large  vortex  structures  travel  downstream  along  the  airfoil  surface.  There  is  no 
vortex  breakdown  observed  in  2D  simulation  since  breakdown  is  a  3-D  and  non-linear 
interaction.  This  shows  2D  direct  simulation  can  only  predict  the  early  instability  stage, 
but  cannot  reveal  the  3D  mechanism  of  flow  transition. 


Figure  13  Contours  of  spanwise  vorticity  from  2D  solution 
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Figure  14  Time  averaged  velocity  field 


Note  that  we  did  not  add  any  disturbance  at  inflow,  there  is  a  question  raised  what 
kind  of  disturbance  triggers  the  instability  wave  inside  the  separated  free  shear  layer? 
Figure  15  shows  time  history  of  fluctuation  pressure  at  different  locations.  We  can  see 
that  the  pressure  fluctuating  appears  first  at  x/c- 1.048  which  is  located  in  the  wake 
very  close  to  the  trailing  edge  (we  call  it  ‘wake  instability’).  As  we  mentioned  above  that 
the  disturbances  introduced  by  initial  field  may  excite  the  most  unstable  instability  wave. 
However,  from  the  simulation  results,  we  find  that  large  pressure  oscillations  appear  first 
in  the  wake  very  close  to  the  trailing  edge.  It  is  because  wake  is  generated  by  a  free  shear 
layer  with  inflection  points  but  without  solid  wall  even  nearby  unlike  the  boundary  layer 
or  separated  mixing  layer,  which  makes  the  wake  most  unstable.  The  time  dependent 
behaviors  of  streamwise  velocity  at  different  streamwise  locations  are  given  in  Figure  16, 
which  also  confirms  that  the  wake  becomes  unstable  first. 


Figure  15.  Time  history  of  pressure  fluctuation  at  different  locations 
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Figure  16.  Time  history  of  streamwise  velocity  at  different  locations 


Then  disturbances  generated  in  wake  near  the  trailing  edge  propagate  upstream 
through  acoustic  waves.  These  upward  traveling  disturbances  could  change  the  pressure 
distribution  on  the  surface.  Then,  the  separation  point  and  the  stagnation  point  become 
unstable  (we  call  it  ‘stagnation  instability’).  The  disturbances  generated  at  these  locations 
converted  downstream  inside  the  separated  shear  layer  as  vortical  disturbances.  The 
shear  layer  has  inflected  profile  with  strong  reverse  flow  which  subjects  to  inviscid  shear 
layer  instability  with  much  larger  growth  rate  than  the  viscous  counterpart.  Figure  17 
shows  the  profiles  for  the  mean  streamwise  velocity  at  different  locations.  Some  DNS 
results  shows  that  the  shear  layer  has  absolute  instability  when  backflows  are  as  high  as 
15-20%  of  external  velocity  (Yang  &  Voke,  2001).  From  Figure  17  we  can  see  that  the 
reverse  velocity  reaches  8%  of  freestream  velocity  at  xIC  =  0.4,  which  is  a  way  below 
the  absolute  instability  range,  but  at  this  point  we  still  can  see  obvious  unsteadiness.  The 
disturbances  are  rapidly  amplified  in  the  shear  layer  and  develop  to  shed  vortex  structures. 
Many  authors  have  related  this  instability  mechanism  to  the  Kelvin-Helmoholtz 
instability.  The  condition  for  the  Kelvin-Helmoholtz  instability  to  occur  is 
0  <  Kh  <  1.2785  (Yang  &  Voke,  2001)  where  K  is  the  wave  number  and  h  is  the  shear 
layer  thickness.  In  the  present  simulation,  unsteadiness  becomes  obvious  at  about 
x/C  =  0.4,  where  the  shear  layer  thickness  is  about  h  ~  0.0202C  and  the  wavelength  is 

X  ~  0.115C .  Then  we  get  Kh  =  — -h  ~  LI 04 ,  which  satisfies  the  criteria  for  the  Kelvin- 

X 

Helmoholtz  instability  to  develop. 
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Figure  17.  Mean  velocity  profiles  at  different  locations 


Figure  18  shows  the  spectrum  of  pressure  waves  at  the  location  close  to  the  wall. 
The  spectrum  at  different  streamwise  locations  from  the  leading  edge  to  the  trailing  edge 
shows  the  strong  peak  at  the  frequency  around  F+  =1.4,  which  is  about  the  vortex 
shedding  frequency  (vortex  shedding  frequency  is  changing  from  time  to  time).  There  is 
a  peak  at  F+  =  0.3  which  may  correspond  to  so-called  low-frequency  flapping.  As  flow 
is  subsonic,  acoustic  waves  travel  in  both  upstream  and  downstream  directions.  The 
pressure  field  is  dominated  by  the  vortex  shedding  frequency,  which  is  also  the  frequency 
of  Kelvin-Helmoholtz  instability  wave.  This  may  provide  a  clue  for  flow  control  that  the 
blowing  jet  for  separation  control  should  have  same  or  similar  frequency  with  very  little 
mass  flow  (sharp  shape  in  both  time  and  space).  When  disturbances  grow  large,  nonlinear 
interactions  take  place.  We  can  see  more  high  frequencies  appear  in  the  spectrum. 
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Figurel8.  Power  spectral  density  of  pressure  at  the  mid-boundary  layer 

4.4.2  Mean  flow 

3D  solutions  are  highly  unsteady.  Mean  flow  characteristics  are  analyzed  first. 
Figure  19  shows  the  maximum  reverse  flow  in  the  wall  normal  direction  along  the 
suction  surface.  The  separated  zone  appears  form  x/C  =  0.19  to  xIC  =  0.68 ,  where  the 
separated  laminar  boundary  evolves  into  reattached  turbulent  boundary  layer.  The  reverse 
flow  reaches  8%  of  free  stream  velocity  at  about  x/C  =  0.5 .  After  that  the  reverse  flow 
appears  to  be  much  stronger.  From  velocity  vectors,  which  are  plotted  in  Figure  20,  the 
separation  zone  can  be  identified  clearly,  which  is  mush  shorter  that  the  2D  sumulation. 
Reverse  flow  is  very  strong  near  the  end  of  the  separation  zone. 
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Figure  19.  Mean  reverse  flow  distribution  on  the  suction  side,  Um,  =  min(w  ) 


Figure  20.  Mean  velocity  vector 

Figure  21  shows  the  mean  pressure  coefficient.  The  flattened  region  indicates  the 
separation  of  the  boundary  layer.  The  strong  adverse  pressure  gradient  at  the  fore  part  of 
the  suction  surface  causes  the  boundary  layer  to  separate  from  the  surface.  In  the 
separation  zone,  adverse  pressure  gradient  is  reduced  dramatically.  Near  the  end  of  the 
separation  zone,  large  adverse  pressure  gradient  forms  again  and  leads  to  a  rapid  increase 
of  reverse  flow  which  is  clearly  demonstrated  in  Figure  19.  Figure  22  shows  the  skin 
friction  coefficient  of  the  suction  surface.  The  skin  friction  increases  when  transition 
happens.  The  separated  shear  layer  becomes  turbulent  and  reattaches  to  the  airfoil  surface, 
forming  a  closed  bubble.  As  the  flow  is  very  unsteady,  the  shape  and  the  length  of  the 
bubble  are  changing  in  a  large  scale.  The  time  averaged  length  of  bubble  is  about 
lb/C  =  0.48  which  can  be  estimated  from  Figure  19.  After  the  reattachment,  the  skin 
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friction  coefficient  grows  rapidly  to  a  high  value  within  a  short  distance,  then  stay  close 
to  that  level  with  some  oscillations  which  are  caused  by  large  vortex  structures. 


Figure  21.  Mean  pressure  coefficient  Figure  22.  Mean  skin  friction  coefficient 

In  this  simulation,  no  external  disturbance  is  introduced,  yet  we  observed  growth  of 
the  disturbance,  instability  waves  and  transition.  We  have  discussed  in  2D  section  that 
according  to  our  observation,  the  acoustic  waves  generated  near  the  trailing  edge  travel 
upstream,  which  could  change  the  pressure  distribution  on  the  airfoil  surface.  The 
stagnation  point  and  the  separation  point  are  closely  related  to  pressure  distribution  and 
can  be  unsteady  if  the  pressure  distribution  changes.  The  disturbances  generated  at  these 
points  are  convected  downstream.  The  separated  shear  layer  is  very  unstable.  The 
disturbances  convected  in  the  shear  layer  are  amplified  rapidly.  The  detached  shear  layer 
becomes  unstable  via  the  Kelvin-Helmholtz  mechanism.  Figure  23  shows  the  maximum 
the  r.m.s.  values  of  fluctuation  velocities  across  the  boundary  layer  along  the  surface.  The 
growth  of  disturbances  can  be  seen  clearly.  Before  the  separation  point  x/C  =  0.19, 
«'and  w'  are  small  but  not  zero,  which  indicates  that  small  2D  disturbances  exist.  The 
disturbances  inside  the  boundary  layer  start  to  growth  before  the  separation  point  (viscous 
instability)  at  around  x/C  =  0.2 .  Then  the  disturbance  grows  in  a  much  faster  rate  in  the 
separated  shear  layer  corresponding  to  the  inviscid  instability  at  x!C  =  0.2  -0.4.  u' 
increases  rapidly  at  about  x/C  =  0.4,  where  the  reverse  flow  reaches  8%  of  the  free 
stream  velocity.  This  strong  reverse  flow  makes  the  shear  layer  very  unstable  and  leads  to 
the  sudden  rapid  growth  of  disturbances.  Following  the  sudden  growth  of  m'.w*  , 
spanwise  velocity  also  starts  to  grow.  The  velocity  fluctuations  reach  a  maximum  value 
before  the  mean  reattachment  point.  The  rapid  growth  of  velocity  fluctuations  indicates 
the  appearance  of  three  dimensional  motion  and  nonlinear  interaction  which  leads  to 
breakdown  to  turbulence.  After  the  breakdown,  the  separated  shear  layer  reattaches  on  to 
the  surface  and  becomes  an  attached  turbulent  boundary  layer. 


38 


0,3 


Figure  23.  Peak  r.m.s  across  the  boundary  layer 

The  mean  velocity  profiles,  the  r.m.s  fluctuation  velocity  profiles,  and  the  Reynolds 
stress  are  depicted  in  Figure  24  (a)  (b)  (c)  respectively.  In  the  separation  zone,  mean 
velocity  profiles  demonstrate  inflected  shapes.  After  the  reattachment,  the  boundary  layer 
develops  to  the  turbulent  velocity  profile.  The  r.m.s  fluctuation  streamwise  velocity 
profiles  at  different  locations  are  displayed  in  (b)  which  clearly  shows  the  maximum 
fluctuation  happened  in  the  separated  shear  layer  away  from  the  solid  surface  in  normal 
direction  corresponding  to  the  invisid  instability  much  larger  than  the  viscous  instability 
near  the  solid  surface.  The  first  location  (x/C  =  0.2)  is  very  close  to  the  separation  point. 
The  peak  appears  at  about  the  center  of  the  boundary  layer.  At  the  following  two 
locations,  the  fluctuations  grow  and  two  peaks  show  up.  From  x/C  =  0.5 ,  the  fluctuations 
grow  rapidly  and  three  peaks  are  found  at  some  locations.  This  evolution  is  related  to  the 
amplification  of  the  upstream  perturbations  due  to  the  existence  of  the  inflected  velocity 
profile  (invisid  instability),  vortex  shedding  and  prime  vortex  breakdown.  Profiles  of 
Reynolds  stress  are  shown  in  (c).  Peak  values  appear  inside  the  boundary  layer. 


Ufljcc 

(a)  Mean  streamwise  velocity  profiles 
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(Maximum  fluctuation  happens  in  the  separated  shear  layer) 


(c)  Reynolds  stress  u'w' 


Figure  24  Statistic  profiles 

4.4.3  Instantaneous  characteristics 

Figure  25  shows  contours  of  the  instantaneous  spanwise  vorticity  in  the  middle  (x,  z) 
plane.  As  3D  simulation  is  started  from  2D  solutions,  the  vortex  evolution  and 
breakdown  can  be  clearly  seen  from  these  time  sequent  pictures.  Before  3D  flow  is  fully 
built  up,  instability  waves  growing  in  the  shear  layer  and  corresponding  to  vortex 
shedding  can  be  clearly  identified.  When  run  time  for  3D  simulation  is  long  enough,  real 
3D  motion  is  fully  developed.  Nonlinear  interactions  of  velocity  fluctuations  become 
very  strong  and  lead  to  rapid  fluctuation  growth.  Compared  with  2D  solution  in  Figure  13, 
organized  vortex  shedding  disappears  and  vortex  breakdown  in  very  short  distance  in  3D 
simulations.  It  is  also  shown  that  2-D  simulation  cannot  catch  the  physics  for  either  the 
flow  separation  and  transition  or  flow  control.  The  2-D  mechanism  is  different  from  the 
3-D  mechanism  and  3-D  simulation  has  to  be  conducted  for  flow  separation,  transition, 
and  control. 
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Figure  25.  Instantaneous  spanwise  vorticity  at  different  time 


The  iso-surfaces  of  instantaneous  vorticity  in  three  directions  are  plotted  in  Figure 
26.  The  transition  process  and  breakdown  of  the  rolling-up  shear  layer  are  clearly 
demonstrated  in  (b).  The  vortex  shed  from  the  separated  shear  layer  are  distorted  while 
traveling  downstream.  The  spanwise  vorticity  iso-surface  becomes  rippled  when  3D 
vortex  appears.  Streamwise  vortice  and  wall-normal  vortices  are  shown  in  (a)  and  (c) 
respectively.  The  interactions  of  3D  structures  cause  spanwise  vorticity  iso-surface  to 
break  into  small  pieces,  indicating  the  breakdown  of  vortices.  The  boundary  layer 
becomes  fully  turbulent  after  reattaches.  During  the  transition  process,  the  instability 
wave  grows  rapidly  and  is  companied  by  the  appearance  of  the  three-dimensional 
motions  which  leads  to  the  breakdown  of  the  boundary  layer. 
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(a)  streamwise  vorticity 


(b)  spanwise  vorticity 
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(c)  wall-normal  vorticity 
Figure  26.  Iso-surface  of  instantaneous  vorticity 


From  the  graphic  (Figure  26)  and  movie  made  with  DNS  data,  we  can  find  the  shed 
vortex  is  quickly  deformed  and  stretched.  The  core  of  the  prime  vortex  is  faded  due  to 
dissipation  and  a  negative  vortex  beside  the  prime  vortex  is  induced  by  the  prime  vortex 
and,  eventually,  the  prime  vortex  breaks  down  to  small  pieced  corresponding  to  transition 
to  turbulence.  From  the  3-D  graphic  (Figure  26),  the  vortex  structure  clearly  shows  the 
stremwise  and  spanwise  wortex  interact  and  a  X  -  shape  vortex  develops,  rolls  up,  and 
breaks  down.  From  Figure  26,  we  can  find  the  streamwise  vortex  plays  a  key  role  in  the 
transition  process.  Whenever  the  streamwise  vortex  appears,  the  vortex  rolls  up  and 
breaks  down  quickly.  There  is  a  remaining  question,  where  is  the  stremwise  vortex  from 
since  there  is  no  3-D  external  perturbation  added?  The  possible  way  could  be  3-D 
acoustic  wave  from  wake  (wake  instability  is  3-D  dimensional)  or  3-D  motion  of  the  shed 
prime  vortex. 

Figure  27  shows  streamwise  evolutions  of  the  vorticity  extrema. (Ox,coy,coz  denote 
vorticity  component  in  streamwise,  spanwise  and  wall-normal  directions  respectively. 
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The  solid  lines  correspond  to  maximum  values,  the  dashed  lines  to  minimum  values. 
The  streamwise  and  wall-normal  vorticity  components  remain  almost  zero  up  to 
transition.  Downstream,  the  minimum  and  maximum  vorticity  oscillates  seriously, 
showing  a  symmetrical  pattern.  The  spanwise  vorticity  component  has  large  value 
close  to  the  leading  edge  and  remains  at  a  lower  level  after  the  separation.  Oscillation 
appears  when  transition  starts. 


x/c 

Figure  27  Instantaneous  extrema  of  vorticity 

Power  spectral  density  of  streamwise  velocity  at  mid-boundary  layer  is  given  in 
Figure  28.  Peaks  at  low  frequencies  located  atF+  =  1  ~  2 ,  as  we  discussed  in  2D  section, 
correspond  to  the  Kelvin-Helmoltz  instability  wave  and  vortex  shedding.  Stronger  peaks 
at  higher  frequencies  appear  gradually  further  downstream  as  a  result  of  nonlinear 
interactions.  During  and  after  the  transition,  energy  peaks  appear  over  a‘  wide  band  of 
frequency. 
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requency  (U  PC) 


The  transition  to  turbulence  on  an  airfoil  with  attack  angle  will  reattach  the 
separated  boundary  layer.  Though  no  external  disturbances  are  introduced,  the  initial 
perturbations  may  come  from  the  up-ward  traveling  acoustic  waves  generated  in  the  wake. 
The  wake,  which  is  a  free  shear  layer  with  invicid  instability,  is  most  unstable  and 
becomes  unstable  first.  The  separation  points  and  the  stagnation  point  are  closely  related 
to  the  pressure  distribution  on  the  airfoil  surface.  The  upward  traveling  acoustic  waves 
may  effect  the  pressure  distribution  which  further  changes  the  location  of  separation 
points  and  stagnation  points.  This  oscillation  introduces  perturbations  inside  the  boundary 
layer.  The  perturbations  are  convected  downstream  as  vortical  disturbances.  The 
separated  shear  layer  has  an  inviscid  instability  and  the  perturbation  will  be  quickly 
amplified  in  a  rate  much  higher  than  the  viscous  instability.  The  traveling  disturbances 
trigger  the  instability  waves  which  are  identified  as  Kelvin-Helmholtz  instability.  The 
appearance  of  3D  motions  of  the  shed  prime  vortex,  where  the  stremwise  vortex  appears, 
and  nonlinear  interactions  of  disturbances  lead  to  the  sudden  growth  of  disturbances  and 
the  generation  of  high  frequencies.  The  breakdown  then  happens  due  to  the  interaction  of 
spanwise  and  streamwise  vortices.  The  shed  prime  vortex  is  stretched  and  deformed,  the 
core  of  the  prime  vortex  is  faded,  a  companion  negative  vortex  is  formed,  and,  eventually, 
the  prime  vortex  breaks  down  to  small  pieces.  The  shear  layer  becomes  turbulence  and 
reattaches  to  the  surface.  This  process  is  then  further  sustained  by  the  global  feed  back 
mechanisms.  This  analysis  leads  to  the  conclusion  that  reattachment  takes  place  right 
after  the  transition,  which  provides  a  clue  for  flow  separation  control.  In  order  to  prevent 
separation  or  reduce  separation  zone,  we  can  use  external  disturbance  to  trigger  the  early 
transition.  When  unsteady  blowing  is  used  in  separation  control,  the  disturbance 
introduced  by  unsteady  blowing  will  excite  the  inherent  local  instability  wave  and  lead  to 
early  transition  to  turbulence  which  will  reduce  the  separation  zone  by  early  reattachment 
of  the  separated  shear  layer  or  late  separation  of  the  boundary  layer.  The  frequency  of  the 
unsteady  blowing  and  the  length  of  the  blowing  hole  provide  the  frequency  and  wave 
length  scale  for  picking  up  the  instability  waves.  Therefore,  the  frequency  of  the  pulsed 
jets  should  coincide  with  the  frequency  of  the  vortex  shedding. 
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Chapter  5 


DNS  for  Flow  Separation  Control  with  Pulsed  Jets 

The  problem  set  up  is  same  as  in  Chapter  4,  but  blowing  jets  have  been  added  in  the 
forehead  of  airfoil  before  the  flow  separation. 

5.1  Results  and  Analysis 

To  study  the  separation  and  transition  processes  on  the  airfoil  and  the  effects  of 
different  types  of  blowing  on  separation  and  transition,  we  set  up  four  cases:  I.  Baseline 
case  without  blowing;  II.  Pulsed  blowing;  III.  Blowing  velocity  with  30°  pitch  angle  and 
90°  screw  angles. 

All  simulations  are  carried  out  with  a  time  step  equal  to  8.35x10 ”5C/C/„.  The 
corresponding  CFL  number  is  around  400. 

5.1.1  Flow  around  the  airfoil  without  blowing  (Baseline  Case) 

Case  I  is  the  baseline  case  without  blowing.  2D  solutions  are  used  as  the  initial  field. 
For  this  case,  the  time  integration  has  reached  t  =  3.474C /Um .  Mean  values  are  obtained 
by  performing  averaging  in  the  spanwise  direction  and  in  time  over  a  period  of  3 C/t/„ . 
The  details  of  the  results  have  been  give  in  Chapter  4 

5.1.2  Flow  around  the  airfoil  with  a  pulsed  blowing  jet 

Based  on  above  observations  and  analysis  for  baseline  case,  we  set  up  two  cases 
(actually  three  cases,  but  the  steady  blowing  has  similar  results  as  the  pulsed  jet  with 
more  mass  flow  and,  therefore,  is  not  reported  here)  to  investigate  how  the  blowing  jet 
(vortex  generation  jet  or  VGJ)  effects  the  transition  process.  All  flow  parameters  and 
geometry  are  the  same  as  in  baseline  case  except  for  that  blowing  is  introduced  on 
suction  side  of  the  airfoil.  The  unsteady  blowing  is  enforced  from  x0  =0.153  to 

x,  =0.175  ,  which  is  before  the  separation  point  xs=0.19.  The  non-dimensional 

frequency  of  blowing  F*  =  FC/U„  is  set  to  be  2  (we  now  believe  1.4  is  a  better  choice), 
where  C  is  the  chord  length  and  F  is  the  frequency.  The  blowing  velocity  has  90°  pitch 
angle  and  0°  screw  -angle.  Pitch  angle  is  defined  as  the  angle  the  jet  makes  with  the  local 
surface  and  skew  angle  is  defined  as  the  angle  of  the  projection  of  the  jet  on  the  surface 
relative  to  the  local  free-stream.  With  this  configuration,  the  blowing  velocity  is  set  in  the 
wall-normal  direction.  The  blowing  velocity  is  given  as  follows, 

w(x,  y,t)  =  Af  (x,  y)exp[-fc(^-l)2] 

where, 

f{x,  y)  =  [0.5  -  0.5  cos{6x  )][o.5  -  0.5cos(^  )J 
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dx  =  2ji— — — ,  9=271^-,  r  =  t- nT,  T  =  UF\  A  =  0.4,  k  =  12 
^“*0  Ly 

Shape  functions  in  space  and  time  which  are  used  to  define  spatial  distribution  and 
temporal  variation  of  blowing  velocity  are  depicted  in  Figure  29.  k  =  12  is  used  in  this 
case.  This  parameter  can  be  used  to  control  blowing  mass. 


Figure  29.  Shape  functions  in  time  and  space 

The  time  integration  for  unsteady  blowing  case  has  reached  t -3.13d  U„.  Time 
averaging  is  performed  over  four  periods  of  blowing.  Mean  velocity  vectors  are  shown  in 
Figure  31a.  It  is  obvious  that  large  separation  zone  which  is  clearly  seen  in  the  baseline 
case  shown  in  Figure  17  is  dramatically  reduced  (almost  removed).  When  we  look  at 
streamwise  mean  velocity  profiles  of  baseline  case  (Figure  17  and  Figure  31a)  and 
blowing  case  which  are  shown  in  Figure  31b,  we  can  find  that  reverse  flow  completely 
disappear  after  x/c  =  0.4  and  boundary  layer  becomes  attached  to  the  surface  afterwards. 


Figure  30.  Mean  velocity  vectors 


47 


i  S 

U/U* 

(a)  base  line  case 


2  3  4  5  6  7 

U/Ucc 

(b)  blowing  case 

Figure  31  Streamwise  mean  velocity  profiles 


The  reduced  separation  zone  can  also  been  seen  from  Figure  32  in  which  the 
maximum  reverse  velocity  is  depicted  against  streamwise  location.  Under  the  effect  of 
blowing,  reverse  flow  happens  in  very  short  distance  from  jc/c  =  0.2  ~  0.3  .  Compared 
with  Figure  19,  the  reverse  flow  quickly  reaches  maximum  value  which  is  about  7%  of 
free  stream  velocity.  While  in  the  baseline  case,  the  reverse  flow  velocity  gradually 
increases  and  reach  maximum  of  about  8%  free  stream  velocity.  This  difference  effects 
the  growth  of  disturbances  which  will  be  discussed  later.  The  distribution  of  mean 
pressure  coefficient  is  plotted  in  Figure  33.  As  the  separation  zone  is  reduced,  the 
flattened  area  no  longer  exists.  Obviously,  this  pressure  coefficient  distribution  will  not 
improve  the  lift  force.  Temporal  variations  of  lift  and  drag  coefficients  which  are 
averaged  over  spanwise  direction  are  shown  in  Figure  34.  For  our  case,  the  reduction  of 
the  separation  zone  decreases  both  drag  and  lift  forces. 
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Figure  32.  Mean  reverse  flow  distribution  Figure  33.  Mean  pressure  coefficient 
on  the  suction  side,  Urev  =  min(«) 


Figure  34.  Temporal  variations  of  lift  and  drag  coefficients 

The  skin  friction  coefficient  distribution  on  the  suction  side  is  shown  in  Figure  35. 
Figure  36  shows  peak  r.m.s  of  velocity  fluctuations.  In  the  separation  zone,  skin  friction 
coefficient  is  negative.  At  around  x/c  =  0A,  there  is  a  rapid  growth  of  skin  friction, 
indicating  the  transition  to  turbulence.  As  discussed  before,  there  are  disturbances  coming 
from  the  leading  edge.  The  rapid  growth  of  streamwise  and  wall  normal  velocity 
fluctuations  starts  after  blowing  injection  point.  After  a  very  short  distance,  they  reach  a 
steady  level.  The  growth  rate  of  spanwise  velocity  is  more  rapid  than  in  baseline  case. 
The  spanwise  velocity  fluctuation  reaches  a  stable  level  around  the  reattachment  point 
x/c  =  0.4 ,  In  this  case,  transition  takes  place  in  very  short  distance  and  disturbance  level 
is  much  lower  than  in  the  baseline  case.  Figure  37  shows  the  r.m.s  streamwise 
fluctuation  velocity  profiles  at  different  x  locations.  Inside  the  separation  zone,  from 
x/c  =  0.2  ~  0.4  ,  peak  values  can  be  seen  inside  the  boundary  layer,  indicating  the 
development  of  the  instability  waves.  After  reattachment,  this  multi-peak  profiles  no 
longer  exist. 
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Figure  35.  Skin  friction  coefficient  distribution  Figure  36.  Peak  r.m.s  across  the 
on  the  suction  side  boundary  layer 


|U*l/Ubc 

Figure  37.  r.m.s  streamwise  fluctuation  velocity  profiles 

The  instantaneous  vorticity  in  span  wise  directions  at  different  time  are  plotted  in 
Figure  38.  The  simulation  is  started  from  2D  solution  as  the  baseline  case.  At  the  early 
stage  of  simulation,  we  can  see  large  vortex  shedding.  The  unsteady  blowing  enforced 
before  the  separation  point  triggers  the  early  transition  of  the  boundary  layer.  The 
boundary  layer  reattaches  shortly  after  separation  and  form  a  short  separation  zone. 
Reattached  boundary  layer  further  becomes  fully  turbulent.  By  comparing  Figure  25  and 
Figure  37,  the  separated  shear  layer  can  be  clearly  seen  in  Figure  24.  In  Figure  38  the 
boundary  layer  is  disturbed  by  the  blowing  and  transition  takes  place  much  earlier  than  in 
the  non-blowing  (baseline)  case  and  the  separation  zone  is  significantly  reduced. 
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t=0.9588 


Figure  41.  Power  Spectral  Density  of  streamwise  velocity 
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From  the  simulation  results  and  analysis  of  this  pulsed  blowing  case,  we  conclude  that 
properly  selected  unsteady  blowing  can  trigger  the  early  transition  through  exciting  most 
unstable  waves  and  non-linear  interactions.  In  this  case,  we  didn’t  observe  the  unstable 
mode  picking  up  mechanism,  because  transition  takes  place  in  very  short  distance.  This 
may  be  due  to  the  large  blowing  mass  flow  which  introduces  a  very  large  disturbance  and 
leads  to  bypass  transition.  To  discover  the  unstable  mode  pick  up  mechanism,  future 
selection  of  blowing  jets  should  be  with  small  blowing  mass  flow  (very  sharp  shape  in 
time  with  large  k)  and  with  a  near  vortex  shedding  frequency.  Though  the  separation  zone 
is  reduced,  both  lift  and  drag  are  decreased.  This  shows  that  the  increase  of  the  ratio  of  lift 
over  drag  can  happen  for  large  attack  angle,  but  not  for  small  attack  angles  (4  for  example 
in  our  case)  although  the  separation  is  almost  completely  removed.  Large  attack  angles  for 
low  Reynolds  number  flow  may  be  selected  as  a  target  for  future  study  of  flow  separation 
control  where  the  ratio  of  lift  over  drag  should  be  significantly  improved. 

5.1.3  The  effect  of  blowing  angle 

To  study  the  effect  of  blowing  angle,  we  setup  another  case  to  simulate  flow  around 
the  airfoil  with  pulsed  blowing  of  30°  pitch  angle  and  90°screw  angle.  The  flow 
parameters  and  geometry  are  the  same  as  the  second  case,  except  the  blowing  velocity 
has  a  pitch  angle  of  30°,  The  vectors  of  blowing  velocity  are  shown  in  Figure  42.  The 
profile  of  the  blowing  velocity  in  streamwise  and  spanwise  directions  are  also  shown  in 
Figure  42.  One  injection  hole  is  centered  at  y/c  =  0.05 . 


Figure  42.  Distribution  of  blowing  velocity 

As  the  time  integration  is  not  long  enough  for  statistic  analysis,  we  only  present 
some  instantaneous  results  here.  More  detailed  analysis  and  comparison  with  the  normal 
blowing  case  will  be  given  in  the  future  paper  or  reports. 
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This  simulation  was  started  from  the  solution  of  the  baseline  case  at  t  =  2.33 C  /  C/„ . 
Figure  43  shows  the  evolution  of  instantaneous  spanwise  vorticity.  At  the  initial  stage, 
the  separation  zone  can  be  clearly  identified.  When  disturbances  introduced  by  the 
blowing  enter  the  separated  shear  layer,  disturbances  are  amplified  when  traveling 
downstream.  At  t  =  0.6253C7 ,  we  can  see  the  rippling  of  the  shear  layer  and  the 
vortex  shedding  after  the  rippling.  At  later  stage,  the  separation  zone  totally  disappeared. 


Figure  43  Instantaneous  spanwise  vorticity  at  different  time 


Figure  44  shows  the  iso-surface  of  instantaneous  vorticity  components  at 
t  =  1.776C/f/„ .  Although  the  experiment  (Bons  et  al,  2001)  shows  the  pitched  jet  will 
obtain  a  better  efficiency  due  to  the  un-symmetric  vortex  structure,  no  much  difference 
from  the  normal  blowing  case  was  found  at  this  simulation  stage.  This  may  be  caused  by 
the  large  blowing  mass  flow  in  both  cases  where  the  unstable  mode  pick  up  mechanism  is 
not  found,  but  bypass  transition  occurs.  The  other  reason  is  the  length  of  integration  in 
time  is  not  long  enough  for  the  statistic  analysis.  More  time  integration  is  needed  for  the 
pitched  jet  case. 
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Figure  44.  Iso-surface  of  vorticity  components 


5.2  Conclusion 

Separation  and  transition  processes  on  a  NACA  0012  airfoil  with  or  without  jet 
blowing  on  the  surface  have  been  investigated.  The  following  conclusions  are  drawn 
based  on  our  observations  and  results  analysis. 


Properly  selected  unsteady  blowing  can  trigger  the  early  transition  through  exciting  most 
unstable  waves  and  non-linear  interactions.  In  our  cases,  we  didn’t  observe  the  unstable 
mode  picking  up  mechanism,  because  transition  takes  place  in  very  short  distance.  This 
may  be  caused  by  the  large  mass  flow  of  the  blowing  jets  which  cause  the  bypass 
transition.  Though  the  separation  zone  is  reduced,  both  lift  and  drag  are  decreased.  Future 
work  should  select  the  large  attack  angle,  small  mass  flow  with  large  k  in  blowing  jets, 
and  appropriate  blowing  frequencies  corresponding  to  vortex  shedding  frequency  which 

is  around  F+  =  1.4  in  our  case. 

The  DNS  data  have  been  post-processed  as  a  movie  which  can  be  obtained  by 
sending  a  request  to  Dr.  Chaoqun  Liu  at  cliu@uta.edu. 
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